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SUMMARY 


Several new methods are developed In this report for 
characterizing atmospheric turbulence, estimating the para- 
meters of these characterizations using modern statistical 
methods, and computing relevant aircraft response statistics 
to such turbulence.* 

In Section 1 of the report, a nongausslan model of atmos- 
pheric turbulence is postulated that accounts for readily 
observable features of turbulence velocity records, their 
autocorrelation functions, and their spectra. New methods 
for computing probability density functions and mean exceedance 
rates of a generic aircraft response variable are developed in 
Sections 2 and 6 respectively. In Section 3y a new method is 
developed for maximum likelihood estimation of tbe parameters 
of a spectrum of known functional form — i.e., the von Karman 
transverse and longitudinal spectral forms. Formulas for the 
variances of such estimates of the integral scale and intensity 
are derived in Section 5. The maximum likelihood method is 
combined with a least-squares approach to yield a method for 
estimating the autocorrelation function parameters of a two 
component model of turbulence in Section 4. Various related 
problems are treated in the Appendices. Implementation of 
the methods to turbulence velocity records is documented in 
another report. 


*Sectlons 1 through M and Appendices A through E of this 
report were completed and submitted to NASA for review 

in May 1978, Section 5 was submitted to NASA in November 
1978 , and Section 6 and Appendices F through K were sub- 
mitted in October I 98 O. 



TURBULENCE MODEL 


In early aeronautical work, turbulence records were 
described as gusts. For purposes of computing wing stresses 
etc., these gusts were modeled In the 1930 's and 19^0 ’s as 
deterministic functions of time — e.g.. Ref. 1. 

During this same period of time, a great deal of funda- 
mental work on the mathematical description of statistical 
phenomena was carried out. Basic mathematical theory of 
stationary time series was developed by Wiener [2,3], Khlnt- 
chlne [4], and Rice [5]. Wang and Uhlenbeck [0] provide an 
excellent description of the state of development In 1945- 
Useful aspects of this work were soon applied to engineering 
problems — e.g., James, Nichols, and Phillips [7]. The 
statistics book by Cramer [5] Is a classic. 

Also, during this same general period, significant pro- 
gress was made by Taylor \_93l0~\ and von Karman [ii] In pro- 
viding a statistical representation of turbulence. Auto- 
correlation functions and power spectra play a fundamental 
role In these representations. An early study of the power 
spectra of turbulence records was carried out by Clementson 
[ 22 ]. 


A principal reason why the power spectrum Is so useful 
a description of random processes Is that It possesses exact 
Input-output relationships for linear time-invariant systems — 
the power spectrum of the output Is the power spectrum of the 
Input multiplied by the square of the magnitude of the system 
frequency response function. Lin [23] generally Is given 
credit for being the first to compute the output (correlation 
function) of a mechanical system from a comparable description 
of Its input. Llepmann [2^] applied these Ideas to the pro- 
blem of buffeting. 

Thus, by the early 1950’ s a methodology existed for 
computing response statistics of an aircraft from a statis- 
tical description of the turbulence excitation provided 
either by the autocorrelation function or Its Fourier trans- 
form, the power spectral density. If the turbulence excita- 
tion Is assumed to be stationary and Gaussian, and the air- 
craft Is modeled as a linear time-invariant system, then 
the response process also Is stationary and Gaussian. In 
this case, the power spectrum of the response provides a 
complete statistical description of the response process. 

However, many turbulence records have a generally non- 
stationary appearance. This nonstationary appearance was 
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taken into account by Press 
records w(t) as homogeneous 
varying standard deviations 


et al. [^5] by modeling turbulence 
Gaussian processes with slowly 
— 1 . e . 5 


w(t) = a(t) z(t) 


( 1 . 1 ) 


where z(t) Is a stationary Gaussian process with zero mean 
value and unit variance, and a(t) Is a nonnegative function 
of time that may be regarded as being either deterministic 
or stochastic. Considerable effort has gone into computation 
of the response statistics of aircraft using the above model. 
Much of this work was done by Press [5 5] and Houbolt 
Rice’s famous formula [5] for the expected number of crossings 
of a stationary random function past a specified threshold 
plays a central role in these studies. Sidwell [57] has 
described Eq. (1.1) as the Press model of turbulence. 

The studies of aircraft response statistics using the 
model of Eq. (1.1) Implicitly assume (1) that fluctuations 
in a(t) occur slowly in comparison with those of z(t), and 
(11) that variations In a(t) are negligible over durations 
comparable with aircraft Impulse-response-function durations. 
Conditions for the validity of assumption (1) are given by 
Mark and Fischer [55] and those for the validity of (11) are 
given by Mark in [55]. When these conditions are satisfied, 
one requires only the power spectral density of z(t) and the 
probability density function of a^(t) to determine the level 
crossing rates and probability density function of an air- 
craft response variable of Interest. The most pertinent 
information about the probability density of a^(t) is con- 
tained In Its mean value and second moment [55]. 

One cannot help but Inquire how far we can progress 
toward computing aircraft response statistics by dropping 
the model of Eq. (1.1) and assuming only that the turbulence 
constitutes an arbitrary (generally nonGaussian) stationary 
random process. This question was addressed by Mazelsky [25] 
in 195^ and slightly earlier In the Russian literature — in 
a more general context — by Kuznetsov and his associates [25]. 
Mazelsky and Kuznetsov showed that higher-order autocorrela- 
tion functions which are time-averaged lagged products ob- 
tained by multiplying turbulence records by themselves 
three, four, and more times also possess exact input-output 
relations — as do their various multidimensional Fourier 
transforms. However, computation times, problems of statis- 
tical confidence of estimates of these characterizations 
obtained from turbulence records of finite duration, and 
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problems of interpretation collectively make these higher- 
order correlation functions and spectra much less attractive 
for the characterization of measured records than the con- 
ventional autocorrelation function and power spectral density. 

Our approach^ therefore is to postulate a turbulenee 
model sufficiently general to include all readily observable 
features of measured turbulence records important to aircraft 
responses ^ and to develop turbulence characterizations from 
this model that have input-output relations sufficiently 
general to predict output probability density functions and 
threshold crossing rates for arbitrary aircraft response 
variables . A requirement of these turbulence characteriza- 
tions absolutely essential for practical application is that 
it be possible to generate realizations of these characteriza- 
tions from measured turbulence records. 


Three Component NonGaussian Turbulence Model 

The simplest model of turbulence as a stochastic process 
is that of a stationary Gaussian process. The vertical 
record shown in Fig. 1 (Ref. 22) Illustrates a turbulence 
record that would appear to be reasonably well modeled as a 
stationary Gaussian time history — especially the portion of 
the record from 120 to 270 sec elapsed time. We might 
reasonably model this record using Eq. ^.1) with a(t) taken 
to be a constant. However, most records of atmospheric tur- 
bulence have a general appearance that is closer to that 
shown In Fig. 2 (Ref. 22). Notice, for example, that the 
portion of the vertical record between 135 and 1^5 sec elapsed 
time has a relatively small rms value; whereas, patches with 
much larger rms values occur shortly thereafter in the 
neighborhoods of 150 and l60 sec elapsed time. Such behavior 
cannot be modeled by a stationary Gaussian process, but can 
be reasonably modeled by Eq. (1.1) when a(t) Is allowed to 
depend on time. 

Each of the records shown In Pig. 2 also exhibits an 
additive weak low-frequency component that appears to fluc- 
tuate independently of the occurrence of the patches. For 
example, during the 5-sec Interval between I 83 and I 88 sec 
elapsed time on the vertical record, high-frequency fluc- 
tuations are absent; however, there remains In that Interval 
a fluctuating weak low-frequency component. Similar but more 
pronounced behavior of this type occurs between approximately 
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FIG. 1. TURBULENCE RECORDS WITH REASONABLY STATIONARY GAUSSIAN BEHAVIOR. 
(Ref. 22, Fig. 17,32. p. 223.) 
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FIG. 2. TURBULENCE RECORDS WITH MILD PATCH-LIKE NONGAUSSIAN BEHAVIOR. 
(Ref. 22, Fig. 17.28, p. 219.) 
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96 and 99 sec elapsed time on the vertical record shown In 
Fig. 3 (Ref. 22). High-frequency fluctuations are absent In 
that Interval; however, there exists there a relatively strong 
low-frequency component. Such behavior cannot be modeled by 
Eq. (1.1). Many other excellent records showing similar 
behavior may be found In Ref. 22. 

The above discussion suggests that we add a low-frequency 
component to the model of Eq. (1.1) — i.e., that we postulate 
a turbulence model [IP] of the form 


w(t) = w^(t) + w„(t) (1.2a) 

5 I 

= ^g(t) + a^(t) z(t) , (l-2b) 

where 

w^(t) = a^(t) z(t) , a^(t) ^ 0 , ( 1 . 3 ) 

and 

E{z(t)} = 0 , E{z^(t)} = 1 . (1.^) 


In Eq. (1.2a), Wg(t) Is the "slow" (low-frequency) component 
and wf(t) is the "fast" intensity modulated component described 
by Eq. (1.3). Since well-behaved turbulence records such as 
the portion of the vertical record shown In Pig. 1 between 
120 to 270 sec elapsed time generally have Gaussian (first- 
order) probability density functions, we shall further assume 
that the stochastic process z(t) In the above model Is sta- 
tionary and Gaussian. In some of the work to follow, we 
also shall assume that V7s(t) Is a stationary Gaussian process. 
Thus, the processes (ws(t)} and {z(t)> are fully described 
by their power spectral densities or autocorrelation func- 
tions. We also shall generally assume that Uf(t) Is a 
stationary random process; however, since af(t)>0, we shall 
not assume that ap(t) Is Gaussian. Furthermore^ the three 
processes {Wg(t)>, {af(t)>, and {z(t)> will be assumed to be 
statistically independent. 

Each of the turbulence records shown in Fig. 4 (Ref. 23) 
clearly illustrates the three turbulence components Wg(t), 
cf(t), and z(t) In the model of Eq. (1.2). Notice that 
throughout each entire record a strong low-frequency component 
Wg(t) Is present. However, for example, between 9 min 0 sec 
and 9 min 45 sec on each record, Wf(t) is negligible In 
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FIG. 3. TURBULENCE RECORDS WITH ADDITIVE LOW-FREQUENCY COMPONENTS PRESENT. 
(Ref. 22, Fig. 17.34, p. 225.) 










comparison with WgCt); whereas, between about 9 min 45 sec 
and 10 min 45 sec, wp(t) grows and then decays back to a 
small value again. Such behavior Is controlled by the temporal 
variations In af(t). None of the three records shown In Fig. 

4 could satisfactorily be modeled by a single stationary 
Gaussian process — or by the process model of Eq. (1.1). 
Equat'lon (1,2) is the simplest model capable of describing 
the overall turbulence behavior illustrated in Fig, 4, 

Locally stationary assumption. For much of the work to 
follow, we shall further assume that fluctuations In Of(t) 
occur slowly In comparison with those of z(t). The quanti- 
tative statement of this assumption [15,15] must be made in 
terms of the power spectral density of the process {z(t)}. 

Let be the two-sided power spectral density of {z(t)}, 

which we define here as the Fourier transform of the auto- 
correlation function of {z(t)} — l.e.. 


^ dx , 

_oo 

where 

(j) (t) = E{z(t) z(t+x)} . 

( 2 ) 

Also let $2 be the second derivative of 


$ 


( 2 ) 

z 


(f ) 


df^ 


^ (f) 
z 


(1.5) 


( 1 . 6 ) 


(1.7) 


Then our locally stationary assumption may be expressed as 


d^£n a^(t) 
dt^ 


<< 32 tt 


^ (f) 
z 


( 1 . 8 a) 


The requirement of Eq . (1.8a) Is derived and discussed In 
Ref. 18 and Is further discussed In Ref. 19. See, In parti- 
cular, Eq. (4.53) and pp . 43 to 53 of Ref. I 8 . When 
has the form of the von Karman transverse spectrum, with 
Integral scale L 25 It is shown in p. 51 of Ref. 19 that the 
requirement of Eq. (1.8a) reduces to 
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(1.8b) 


d^iln a^(t) 
dt^ 


< 0.08 — 

2 


T 


3 


where V is the aircraft speed. This requirement is more 
meaningfully expressed in terms of written as a function 
of the spatial variable x = Vt : 


d^£n a^(x/V) 
dx^ 


0. 08 



(1.8c) 


The examples of this requir>ement discussed on pp . 43 to 30 of 

Ref. 18 show that when a typical length sea le of the func- 
tion Oj^(x/V) is at least 10 times the locally stationavy 

assumption of Eq. (1.8) is well satisfied. 

A more relaxed locally stationary condition expressed In 
terms of the autocorrelation function of £n a|>(t) and the 
aircraft frequency response function of interest is given by 
Eq. (5*17) on p. 56 of Ref. 19. An expression for the auto- 
correlation function of in a|‘(t) in terms of measurable 
quantities is given by Eq. (5.36) on p. 6l of Ref. 19. 

For typical atmospheric turbulence records, the locally 
stationary assumption of Eq. (1.8) is believed to introduce 
negligible error In the results to follow. The simulation 
studies carried out In Ref. 2h support this conclusion. 

Spectral form of z(t). In some of the work to follow, 
we shall assume that the process {z(t)} has the appropriate 
(transverse or longitudinal) von Karman spectral form. When 
the locally stationary condition of Eq. (1.8) Is satisfied, 
it is shown In Ref. l8 that the spectral form of wf(t) = Of(t)x 
z(t) Is unaffected by the fluctuations in Cf(t) — l.e., 
wf(t) will have the same form of spectrum as {z(t)}. 

Comparison with previous related models. A three- 
component turbulence model functionally similar to Eq. (1.2) 
has been studied by Reeves et al . [25—27] and Sldwell [25]. 

However, both Reeves and Sldwell assume that their counter- 
part to our Qf(t) is a Gaussian process with zero mean value. 
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This restriction forces their counterpart to Wf(t) to have 
periods of very low Intensity [when cTf(t) Is near its mean 
value of zero]. Reeves motivation [25—27] for adding his 
counterpart to wg(t) is to partially remove such ”deep fades.” 
Thus 5 in Reeves' model, the power spectra of his counterparts 
to Ws(t) and Wf(t) are taken to have the same Dryden form. 

On the other hand, our approach in Refs. l8 and 19, and 
in the present work, is to introduce a minimum of assumptions 
pertaining to the behavior of Cf(t) and Wg(t) and to extract 
descriptions of these processes relevant to the aircraft 
response problems from measured turbulence records. Algo- 
rithms for computing the power spectra of z(t), Ws(t), and 
a|>(t), and for computing moments and probability density 
functions of Ws(t) and af(t) are provided in Sec. 6 of Ref. 

19. Some of these techniques are modified and extended in 
the present work. 


Aircraft Response Metrics 

Aircraft model. It is our goal to develop expressions 
for the (first-order) probability density functions and the 
threshold mean crossing rates for a general aircraft response 
variable — using the above described model of turbulence as 
the excitation or input. We shall model the aircraft as a 
linear two-terminal time-invariant system described either 
by its unit-impulse response function h(t) or complex frequency 
response H ( f ) : 

h(t) dt . (1.9) 


For any turbulence sample function w(t), the aircraft response 
y(t) is the convolution of w(t) and h(t) — i.e., 


y (t ) 


0 00 

h(x) w(t-x) dx 

— 00 


.oo 

w(x) h(t-x) dx 

— 00 


(1.10a) 


(1.10b) 
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Implicit In this treatment Is the assumption that the spatial 
variation of w(t) Is negligible over dimensions comparable to 
those of the aircraft. The aircraft Impulse response may 
represent the response of any aircraft variable of Interest — 
e.g., the stress at a critical point in a wing span — to the 
turbulence velocity Input w(t). The impulse response h(t) 
may Include the action of a pilot or autopilot modeled as a 
linear feedback element. That Is, h(t) may be thought of 
as either an open or closed loop unit-impulse response func- 
tion . 


AutoooTvetation function input-output relationship . It 
Is well known — e.g., p. 71 of Ref. 29 — that the power 
spectral density ty(f) of the aircraft response y(t) Is re- 
lated to the Input power spectrum ^-^(f) and aircraft frequency 
response H(f) by 

$„(f) = |H(f) |2 . (1.11) 

y w 


Let c1)^(t) and c{)y(T) denote, respectively, the autocorrelation 
functions of the excitation and response processes {w(t)} 
and {y(t)}: 


9 ( ) 


<fy(T) 


$ (f) e df 

J w 

— 00 
^00 

$ (f) 

J y 

— 00 


( 1 . 12 ) 

(1.13) 


Then, from the convolution theorem and Eq, (1.11), it follows 
that the autocorrelat Ion function Input-response relation- 
ship is the convolution 


♦ y(T) 



— 00 


(t-O 

(C-t) 
(C+t ) 






(l.l4a) 

(1. l4b) 

(1. l^c) 
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where 


= J 


H(f) 


12Trf T 
e 


df , 


( 1 . 15 a) 


0OO 

h(t) h(t+T) dt , (1.15b) 

J 

— 00 


and where Eq . (I,l4b) follows from Eq. (I.l4a) and the fact 

that the autocorrelation function 4)h(T) Is an even function 
of T, and Eq. (I.l4c) follows from the fact that (|)y(T) Is an 
even function of t. 


Since we have assumed that the turbulence components 
{Wg(t)}, {af(t)>, and {z(t)} are mutually statistically In- 
dependent ^ It follows that {ws(t)> and {wf(t)}also are mutually 
Independent (and therefore uncorrelated). Therefore, from Eq. 
(1.2) It follows that 


(p (t ) = (^ (t ) -h (p (t) 
s f 


(l.l6a) 


= (t) + (p (t) (P^(t) 

s f 


( 1 . 16 b) 


where the second line Is a consequence of Eq. (1.3) and the 
assumed Independence of {a^(t)l and {z(t)>. Moreover, from 
the locally stationary assumption of Eq . (1.8) It follows 

\_ 18 ~\ that 


= E{ap ; 

hence, we have from Eqs. (I.l6c) and (1.17c) 


( 1 . 17 a) 

(1.17b) 

( 1 . 170 ) 

( 1 . 18 ) 
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That is^ our turbulenoe model of Sea. 1.1 implies that the 
autocorrelation function of the turbulence {w(t}} will appear 
as a superposition of the autocorrelation function of the 
slow component Wg(t) and an amplitude scaled version of the 
autocorrelation ^^(t)of the Gaussian component z(t). 


The behavior of our turbulence model described by Eq. 

(I.l8) Is Illustrated in Pig. 5 for idealized autocorrelation 

functions. Notice that (t) Is shown decaying much more 

s 

slowly than (|) 2 (t). This latter behavior is consistent with 
our assumption that the process Wg(t) fluctuates slowly in 
comparison with the process z(t) — as Is evident from the 
record shown in Pig. 4. In fact, Wg(t) may be thought of as 
a slowly varying mean wind; whereas, wp(t) = Of{t) z(t) may 
be thought of as ordinary turbulence. 


Substituting Eqs. (I.l6 to l.l8) into Eq . (I.l4c), we may 
express the aircraft response autocorrelation function as 


( T ) 

y 


— C» 
^00 

-GO 

C 00 


[^w ^w 

(C) + E{ap 
s 


( 1 . 19 a) 

(1.19b) 

(1.19c) 


which are the desired autocorrelation function input-response 
relationships. Equations (1.19a) and (1.19b) are an exact 
consequence of the turbulence model described by Eq. (1.2) 
and the assumed Independence of {wg(t)}, {cf(t)}, and {z(t)}; 
whereas. In Eq. (1.19c), the locally stationary assumption 
of Eq. (1.8) has been used. 

Aircraft mean-^square displacement and velocity responses . 
Setting T = 0 In Eq. (1.19) directly yields the mean-square 
aircraft displacement response in terms of the turbulence 
component autocorrelation functions: 
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E{yM 5 


$y(0) 


.00 


j 

— CO 


C'J’w <5.5 


r°o 


4>v,(€) [<(>„ (5) + (5) ip^U)2 d5 


a. 


J [cj)^ (^) + E{ap (f)^(^)] , 


w 


'‘'z 


(1.20a) 

(1.20b) 

(1.20c) 


where the approximation of Eq. (1.20c) again depends on the 
locally stationary assumption of Eq. (1.8). 

Comparable expressions may be written for the aircraft 
mean-square velocity response E{y^}. If we formally dif- 
ferentiate Eq. (1.10), we obtain 


y(t) 


^00 

w(t) h(t-x) dT , 

J 

~co 


( 1 . 21 ) 


where h(t) Is the time derivative of the displacement Impulse 
response h(t). Since some h(t) of Interest may contain dis- 
continuities — e.g., at t = 0 — care must be taken In computing 
i^(t). That is, fi(t) must satisfy 


h(t ) 



— 00 


h(^) d^ 


( 1 . 22 ) 


Thus, If h(t) has a discontinuity then fi(t) must contain a 
delta function at the same place so that Eq, (1.22) is 
satisfied . 

As In Eq. (1.15b), we may define for h(t) 

I CO 

h(t) h(t+x) dt . ( 1 . 23 ) 
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The mean-square velocity response may then be expressed in 
the same manner as in Eq. (1.20): 


E{y"} E c|).(o) 



1 — 1 
-e- 
s: 

+ 

w. 


S 


h 



+ 

4> 


s 


f 


(C) 

+ 

E{a 


^)] (1.24a) 

C) dC (1.24b) 

} , (1.24c) 


where the approximation of Eq. (1.24c) again depends on the 
locally stationary assumption of Eq. (1.8). 


Parametric Description of Atmospheric Turbulence 

In order to tabulate relevant features of measured tur- 
bulence records. It is necessary to characterize such records 
by a set of parameters that can be extracted from the records. 
The Integral scale and mean-square value of a record are 
examples of such parameters. 

ChavacteT'lzat'lon of "stow'^ component Ws(t). The NASA MAT 
program (Measurement of Atmospheric Turbulence) \_Z0~\ has 
concentrated on obtaining atmospheric turbulence recordings 
accurate to frequencies (wavenumbers) well below typical 
positions of the "knees" of von Karraan spectra. This effort 
has required exceptional care In aircraft Instrumentation [3i]. 
Typical autocorrelation functions computed from MAT records 
suggest that an efficient characterization of autocorrelation 
function of fhe low-frequency component Wg(t) Is a low- 

order polynomial approximation to 4>w (?) valid In the neigh- 

s 

borhood of ^ = 0 [pp. 64,65 of Ref. 19]* Although this low- 
order polynomial representation may be Interpreted as the 

first few terms of a Maclaurln expansion of c|)^ (?)a we shall 

s 

use a (constrained) least-squares procedure to compute the 
actual expansion coefficients as described In Sec. 4 of 
this report. 


18 



Figure 6 displays the autocorrelation function of the 
vertical record shown In Pig. According to Fig. 6, about 

JS% of the ’’power" In the vertical record of Fig. 4 Is In 
Ws(t) and about 22% Is In wf(t). 

To motivate our concentration on region of ^ 

near the origin, we note first that to predict the first-order 
probability density functions and threshold mean crossing 
rates of aircraft responses, we require E{y^} and E{y^} [pp. 
34—50 of Ref. 19]- Expressions for these quantities are 
given by Eqs. (1.20)^and (1.24) of this report. If we de- 
compose E{y^} and E{y^} into contributions from our ’'slow’’ 
and "fast" components ws(t) and Wf (t ), respectively , 


E{y^} = E{yg} + E{yp 

and 

E{p} = E{y|} + E{yp , 

then we see from Eqs. (1.20) 
contributions are given by 


(1.25) 

( 1 . 26 ) 

and (1.24) that the slow component 


E{y|> 

and 

E{yU 


j s 

r ^ 

<t>W 

j s 






(1.27) 


(1.28) 


Let us now assume — for the purposes of the present discus- 
sion — that h(t) and h(t) are of duration sec only — e.g., 
that h(t) and h(t) are zero outside the Interval 0;<t<T]^. 

Then, it is easy to show from Eqs. (1.15b) and (1.23y, re- 
spectively that c{)]^(5) and are zero for 1^|>T^. Con- 

sequently, from Eqs. (1.27; and (1.28), we see that E{y|} 

and E{^|} depend on the values of 4)^^ (O only for values 

, s ^ 

of ^ satisfying 
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VERTICAL GUST VELOCITY 
AUTOCORRELATION 



0123456789 xIO* 

f (meters) 


FIG. 6. 


AUTOCORRELATION FUNCTION OF VERTICAL RECORD SHOWN IN FIG, 4. [MOUNTAIN 
WAVE CONDITIONS. AIRCRAFT SPEED 197 m/sec (646 ft/sec.)] 



As an example [2J], we consider an aircraft flying at 
Mach 2.7 at an altitude of l 8 km (60,000 ft) where the speed 
of sound Is approximately 297 m/sec (975 ft/sec). If we take 
Tk = 3 sec, then the correlation interval of Interest Is about 
2400 m (7900 ft). On the other hand, if we take = 10 sec, 
the correlation Interval of Interest is about 8000 m (26,200 
ftK For the autocorrelation function displayed In Pig. 6 , 
we see that, for the 3-sec Impulse response, a linear approxi- 
mation to (j)iw provides a good fit; whereas, for the 10 -sec 
s 

Impulse response, a quadratic approximation would be adequate. 
Consequently , we shall take for our parametric representation 

*|>w 4*w 

^ s s 

m 

"fw (5) = I a s-’ , (1.29) 

~"s J=0 J 


where the degree m chosen for the above polynomial representa- 
tion of (^) is to depend on the observed complexity of 
s 

(5) and the Interval in ^ over which (O is to be 
s s 

represented by Eq. (1.29). 

Two fundamentally different approaches may be used In 
generating the representation of Eq. (1.29). On the one hand, 
we may take the expansion Interval as 0<_^<_T]^ and Include both 
odd and even powers of ^ in the expansion. This procedure 
clearly is the best for the autocorrelation function shown in 
Fig. 6 . On the other hand, we may take for the expansion 
interval, the even Interval -I'h^^l.l'h- Since, by definition, 

(})w ( 5 ) must be an even function of 65 this latter approach 
s 

must contain only even powers of Generally, the latter 

approach will require higher powers of 4 to get a good fit 
using Eq. (1.29) — l.e., a larger value of m — but it has 
the advantage that the integrals obtained by substituting 
Eq. ( 1 . 29 ) into Eqs. (1.27) and (1.28) are. In some situa- 
tions, more easily evaluated. 

The reader may wish to Interpret Eq. (1.29) as a trun- 
cated Maclaurln series expansion. In the first of the above 
two approaches where Eq. (1.29) applies to the Interval 

0 <_ 5 £Thj the derivatives of (^) in the interpretation must 

s 

be considered as one-sided derivatives valid only In the 
region ^^0 — l.e., we have aj = cj) ^^* ^ ( 0 + )/j ! , where c))^^*^( 0 +) 
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denotes the jth order derivative evaluated at ^ = 0+. In 
the second of the above two approaches where Eq. (1.29) ap- 
plies to the interval we have aj = ^ ( 0 )/j ! , 

where here no "one-sided condition" is required since ^ in this 

interpretation^ it is assumed that is continuous at 

^ = 0 for j = 0,1,2,. ..^m. 


Since the integrands of Eqs. (1.27) and (1.28) are 
necessarily even, our estimates of E{y|} and E{;ys} obtained 
using Eq. (1.29) may be expressed as 


and 


E{y"} == 2 I 

® j=o 





E{y|} 


m 

I 

J=0 


a 


J J 




dS 


(1.30) 

(1.31) 


Equations (1,30) and (1.31) are valid for either of the above 
types of expansion. However, from the Fourier mate to Eq . 
(1.15a), we have 


• oo 

H(f) = J ,fj^(5) d5 


Differentiating Eq. (1.32) j times, we find 


( 1 . 32 ) 


^fJ J ^ 

.00 


hence, setting f = 0 in Eq. (1.33), we have 


■” 5 J $.(€) dS = — ^ 

(-12tt)'^ df'^ 


J 

.00 


(1.34) 


f=0 


For J = odd, both sides of Eq. (1.34) vanish; however, for 
j = even, we have 


( j )^(0 d ? = 2 


^ J = even. 


( 1 . 35 ) 
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Thus, the terms j = even in Eq. (1.30) may be evaluated from 
|H(f)|^ using Eqs. (1.34) and (1.35). Furthermore, for well- 
behaved h(t), we have from the Fourier mate to Eq. (1.9): 


h(t) 


• CO . 

i27rfH(f) df ; 

J 

^OO 


( 1 . 36 ) 


hence, i2TrfH(f) is the Fourier transform of h(t). It follows 
from this fact and Eq. (1.34) that 


cj)^(^) d^ 


47t^ d’^' 

(-12tt)^* df'^' 


fH(f ) 


f=0 


(1.37) 


For j = odd, both sides of Eq. (1.37) vanish; however, for 
j = even, we have 

-oo ^ ^oo 

(f)^(?) dC = 2 1 d^ , j = even. (1.38) 

_oo 0 


The terms J = even in Eq. (1.31) may be evaluated from |H(f)|^ 
using Eqs. (1.37) and (1.38). 

For situations where odd powers of j are Included in 
Eqs. (1.30) and (1.31), a different approach is available 
for evaluating these expressions. We may decompose cf)y(T) 
into components arising from Ws(t) and Wf(t) where, from 
Eq. (1.19a), we see that the contribution from Wg(t) can be 
written as 


S 


j (j)j^(C+T) (j)^ (K) d^ 
s 


(1.39) 


Dif f erentlating Eq . (1.39) twice yields 


(t) = 
^ s 


(^) dC 

s 


(1.40) 
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where we note that 


E{yg} = ( 0 ) = -j 

3 


(5) <55 , 

s 


(1.^1) 


see - e.g., p. 33 of Ref. l 8 . In regard to the behavior of 

at ^ = 0 , we point out that inclusion of mass or inertia 
in the aircraft impulse response h(t) guarantees continuity of 
at C = 0. Since 0^(5) is an even function of 
also must be even; hence, we may rewrite Eq. (l.Al) as 


E{y 


= -2 j 4,^(5) 4„ (5) dC , 




While, from Eq. (1.39):, we may express E{y|} =: (0) as 


E{yg} = 2 ^w 


(1.^3) 


Combining Eqs . (1.43) and (1.42) with Eq. (1.29) yields 


m /•«> . 

E{yU ^ 2 I 5'^ 4v,(5) ds , 

s J=o J i 


(1.44) 


and 


E(yp 


-2 


m 

I 

j=0 


5 J 4" (5) d5 , 


(1.45) 


where Eqs. (1.44) and (1.45) are valid for cases where both 
odd and even powers are included in Eq. (1.29) 5 or cases 
where only even powers are included. Furthermore, it is 
possible to avoid taking moments of the derivatives of (pj^(^) 
as we now show. Using integration by parts, we have 


5 '’ 4{;(e) d5 = S'’ 4^(S) 


-J 


/■«> . 

S'’ 5 4,-(5) ds 


S^"’'- 4^(0 dS , j i 1 , 


(1.46) 
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since the first term in the right-hand side of the first line 
vanishes for J^l. Repeating the integration by parts on the 
right-hand side of Eq. (1.46), we further have 




dK = -3 


-j-1 




-(J-1) 


5^ 




= J(J-l) 


■J-2 


<l>h(5) 


d? 


J > 2 




since the first term in the right-hand side of the first line 
vanishes in this case for j^2. Combining Eqs . (1.46) and 

(1.47) with Eq. (1.45) yields 


E{y^} 


-2 


^OO rCO 

<(>f;(5) d5 - a, 1 d5 

0 0 


+ 'I J (j-1) a. 
J = 2 ^ 


5^ ^ <).j^(5) d5 . 

0 


(1.48) 


Equations (1.44) and (1.48) are particularly suitable for 
computing E{y|} and E{y|} in cases where the "one-sided" 
expansion Including odd powers of ^ is used in Eq. (1.29). 

In view of the fact that the first or higher order derivatives 
of turbulence autocorrelation functions may not be continuous 

at ^ = 0, representation of (5) over the Interval 

s 

by Eq. (1.29) using odd as well as even powers of and then 
computing E{yg} and E{y|} with Eqs. (1.44) and (1.48) is 
probably the best overall method. It is easy to relate the 
moments in the right-hand sides of Eqs. (1.44) and (1.48) to 
derlvates of the unilateral Laplace transform of cj)j-^(5) — 
which may be useful in evaluating the moments. 

The autoaorretat'Con funct'ion (K) of the ''stow'* oom- 

s 

ponezit WQ(t) oonta-ins complete Information about the power 
spectrum of WQ(t)^ because the two are a Fourier transform 
pair. The most useful general set of parametric descriptors 
of this Fourier transform pair appears to be the set of co- 
efficients ajj g=0^1y...^m of a power series representation 

of ^w (K)' Blnce (^) Is necessarily an even function of 
s s 
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th'is r*epv esentati-on may be chosen to appToxtmate (^) ovev 

s 

etthev or where Tj^ is the duration of the 

aircraft impulse response of interest. The representation 
over 0<^^<fFf^ generally will he preferable because of its more 
rapid convergence. In this case^ Eqs. (1.44) and (1.48) may 
be used to compute E{y^} and However ^ the representa- 
tion over the interval -Iji—^—^h evaluation of all 

integrals in Eqs. (1.30) and (1.31) directly from the squared 
aircraft frequency response magnitude \H(f)\'^ using Eqs. 
(1.34)^ (1.35), (1.3?), and (1.38). 


In addition to a representation of shall re- 

s 

quire the first few moments of the probability density of 
Ws(t), say the first four, in order to verify the Gaussian 
property of Ws(t). A method for computing these moments and 
forming an approximation to the first-order probability den- 
sity of Ws(t) is described in Sec. 6.4 of Ref. 19. 


Characterization of stochastic intensity Oj-(t) of the 
"fast” component Wf(t). The locally stationary assumption 
of Eq. (1.8) of this report Involves only turbulence charac- 
teristics. For this assumption to be satisfied, variations 
in the intensity cTf(t) of the fast component must be small 
over Intervals comparable with the Integral scale or nominal 
correlation interval of the component z(t). Two additional 
locally stationary conditions were studied in Ref. 19* These 
additional two conditions are given by Eqs. (3.43) and (3.46) 
on p. 32 or Eqs. (5.9) and (5.10) on p. 54 of Ref. 19. 
Equations (5.9) and (5.10) are a statement of these conditions 
for engineering purposes, and are somewhat more easily satis- 
fied than Eqs. (3.43) and (3.46). The physical interpretation 
of the combination of these two additional conditions is 
that variations in a^(t) must be negligible over durations 
comparable to the duration Tj^ of the aircraft impulse response 
function of interest . Thus, satisfaction of these additional 
two locally stationary conditions depends on aircraft charac- 
teristics as well as the behavior of the turbulence component 

Uf (t ) . 

For situations where these additional two locally sta- 
tionary conditions are satisfied, it was shown in Secs. 4.1 
to 4.4 of Ref. 19 that the first-order probability density 
function and threshold mean crossing rates of an aircraft 
response variable can be computed from the probability density 
function of and furthermore, that the first few moments 

of the probability density of provide the most Important 
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parametric descriptions of a|.(t). Although a careful study 
of these additional locally stationary conditions Involving 
recorded turbulence data and actual aircraft characteristics 
has not been made, we generally would expect these additional 
two locally stationary conditions to be satisfied for sub- 
sonic aircraft for engineering purposes. However, these 
additional two conditions may not be satisfied in the case 
of supersonic aircraft. 

It will be shown in later sections of this report that 
when the additional two locally stationary conditions are not 
satisfied, we shall require the power spectral density 

f 

or autocorrelation function cr^(t) to predict 

the most Important nonGaussian correction terms for the first- 
order probability density of an aircraft response variable. 

The most important parametric descriptors of are the 

f 

first few "one-sided" or "two-sided" power series expansion 
coefficients of (|)q2(t) — as was the case for c})^ (t). 


Characterization of stationary Gaussian component z(t) 
of the ''fast” component Wj^(t). Since z(t) is, by hypothesis, 
a stationary Gaussian process with zero mean value and unit 
variance, it is completely described by its power spectral 
density or autocorrelation function. In our computational 
work, we shall assume that z(t) possesses the appropriate 
(transverse or longitudinal) von Karman spectral form. For 
these cases, z(t) is completely described by a single para- 
meter — the integral scale of the appropriate transverse 
or longitudinal von Karman spectral form. 


Summary of Turbulence Model Charac ter i za ti ons 

Basic Model 

w(t) = WgCt) + c^(t) z(t) 

a^(t) > 0 , E{z} = 0 , E{z^} = 1 , 

Wg(t) stationary and Gaussian , 

a^(t) stationary. 
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z(t) stationary and Gaussian 
w (t), a^(t), and z(t) mutually Independent. 

o X 


"Complete" Charaeterizat-ion of Model 


Component 


Characteri zat i on 


Wg(t ) 


a^(t ) 

z (t ) 


P.S.D. or Autocorrelation func- 
tion of WsCt), (Also P-D.F. of 
ws(t) to ’’check" Gaussian 
Assumption . ) 

P.S.D. or Autocorrelation Func- 
tion of a^(t), P.D.F. of a^(t). 

P.S.D. or Autocorrelation Func- 
tion of z(t). 


Parametr-ic Chavaoterlzat'ion of Model 


Component 

Wg(t ) 


a^(t ) 


z (t ) 


Characteri zati on 

"One-sided", or "two-sided" power 
series expansion coefficients 
of autocorrelation function of 
ws(t), (Also first few moments 
of P.D.F. of Wg(t) to check 
Gaussian assumption. ) 

"One-sided" or "two-sided" power 
series expansion coefficients 
of autocorrelation function of 
ai(t). First few moments of 
P.D.F. of a^(t). 

Integral scale of appropriate 
(transverse or longitudinal) 
von Karman spectrum. 
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AIRCRAFT RESPONSE PROBABILITY DENSITY FUNCTIONS 


In Sec. 4.4 of Ref. 19a a series expansion was developed 
for the first-order probability density function of an arbi- 
trary aircraft response variable. The expansion Is valid for 
situations where the three locally stationary conditions 
described above and in Ref. 19 are valid. In this section, 
a similar series expansion for the first-order probability 
density of an arbitrary aircraft response variable is devel- 
oped; however, in this new treatment only the first locally 
stationary condition described by Eq. (1.8) is required. This 
condition depends on turbulence properties only, and is 
believed to be virtually always satisfied. Thus, the present 
results apply to supersonic aircraft with arbitrarily high 
Mach numbers — as well as to subsonic aircraft for which the 
simpler results of Sec. 4.4 of Ref. 19 apply. 

Gaussian Property of Response Process Conditioned 
on the Intensity Modulation Process a^(t)* 

In Sec. 2, we shall assume that the "slow" component 
Wg(t) and the component z(t) in the turbulence model of Eq. 
(1.2) both are stationary Gaussian processes with zero mean 
values. When Ws(t) and z(t) satisfy this zero mean and 
Gaussian assumption^ the response process y(t)^ conditioned 
on the process o^(t) is a zero mean strictly Gaussian (gen- 
eraity nonstationary ) process . To prove this, we note first 
that each sample function of the response process (y(t)} can 
be expressed as 

y(t)[a^(u) = YgCt) + y£.(t)la^(u) , -°o<u£t (2.1) 


where yg(t) is the aircraft response to the "slow" turbulence 
sample function w (t) In Eq. (1.2), 

,00 

y^(t) = 1 h(x) Wg(t-T) dx , (2.2) 

0 


*This section closely parallels Sec. 4.1 of Ref. 19* 


29 


and yf(t)|af(u) is the aircraft response to the "fast" tur- 
bulence sample function w^-Ct) In Eq . (1,2), 


(t ) I a^(u) 


,00 

h(T) [w^(t-T ) I a^(u) ] dT, -«3<u<_t 

0 


.oo 

h(T) a-:,(t-T) z(t-T) dx . 

J i 

0 


(2.3a) 

(2.3b) 


The vertical bars followed by a^(u) In Eqs. (2.1) and (2.3) 
denote that. In the stochastic process Interpretation of these 
equations, the sample function a^(u) Is assumed to be known 
or specified for all values of u within the Interval -«><u£t . 
Also, In writing Eqs. (2.2) and (2.3), we have assumed that 
h(t) Is deterministic and causal — i.e., that h(t) - 0 for 
t<0. When ap(u) is assumed to be known or specified for all 
~co<u£t, the right-hand side of Eq. (2.3b) represents a (deter- 
ministic) linear transformation of the Gaussian random function 
z(t). Thus, the conditional random process {yp ( t ) | 0 £*(u) } , 
where Of(u) is specified for all -«><u:<t, is itself strictly 
Gaussian (and generally nonstationary) — since any linear 
transformation of a Gaussian process is itself Gaussian — 
e.g., Cramer [5], pp. 312 and 313. Furthermore, {wg(t)} is 
assumed to be Gaussian; thus, from the linearity of Eq. (2.2), 
the random process {yg(t)} also is Gaussian. Moreover, since 
{wg(t)} and {z(t)} are assumed to be independent, {ys(b)} 
and {yp (t ) I Of (u) } also are independent. Finally, since the 
sum of any number of independent Gaussian processes is neces- 
sarily Gaussian — e.g., Cramer [5], p. 316 — it follows from 
Eq. (2.1) that the conditional response process {y ( t ) | ap(u) } , 
where ap(u) is specified for all -°o<u£t , is strictly Gaussian 
(and generally nonstationary) . This result does not depend 
on any locally stationary assumption. From the zero mean 
value assumptions for the processes {wg(t)} and {z,(t)} it 
follows further from Eqs. (2.2) and (2.3b) that {y(t)|ap(u)} 
also has zero mean value. 

Let us denote the conditional mean square value of the 
process {y(t)} by 

. of = a^(t) = E{y^ (t ) I a „(u) } , -<»<u<t ; (2.4) 

y y ' 1 - ^ 


30 



I 

i 

I 

I 

I 


that Is 5 cTy is the expected value of the squared system re- 
sponse given that af(u) is specified for all -«><u<_t . The 
expectation operation therefore takes place over the ensembles 
of input functions {Wg(t)} and {z(t)} with af(u) being con- 
sidered as known. Since {y(t)|cf(u)} is generally non- 
stationary, Oy is generally a function of t. Furthermore, 
when we (later) consider Of(t) to be a stochastic function 
of time, CyCt) also becomes a stochastic function of time. 

Let us now consider the function Cf (u) , -“<u£t to be 
the limiting case of an infinite dimensional "vector" 

E d£.(ui ,U 2 , . . . ,Uj^) as n is taken to approach infinity 
and uj+i - Uj is shrunk to zero for all j = 1,2,..., n-1. 

Thus, the assumption that the "vector" g^ is specified is 
identical to an assumption that the function af(u) is speci- 
fied for all -°o<u<^t . Let p(y|uf’) denote the conditional 
probability densit’y of the aircraft response y(t) given that 
Cf(u) is specified for all -~<u£t . Then, from the above 
discussion, p(ylqf) is strictly normally distributed with 
variance Cy = CyCt) described by Eq. (2.4): 


n ^ 7 

p(y|g.) = , e ^ . (2.5) 

/2tt 

y 


Series Expansion of Response Probability 
Density Functions 

Equation (2.5) expresses the conditional probability 
density of the aircraft response y(t) given that the random 
function 0f(u) is specified for all values of 0<u£t . The 
unconditional probability density of the aircraft response 
is the expectation of p(y|of) with respect to the joint 
probability density of the'^vector g^.; i.e.. 


P(y) 


#oo 

p(ylSf) p(?f) 

■) 


(2.6a) 


= E{p (y I g^.) } 


(2.6b) 
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where the Integral In Eq. (2.6) represents the limiting case 
of an infinite order Integral over the n-dimenslonal space 
{Uj 3 U 2 , . . . as n is taken to approach infinity and 
^j +1 ” shrunk to zero for all J = 1 , 2 ,... ^n-l as 

is described above. (It will become evident shortly that no 
such Integrals will have to be carried out to apply the 
methods being presented.) 

In most turbulence records , fZuctuat-ions in the stochastic 
function Qf(t) are of the order of not more than 50 % of the 
mean value of Cf(t). Each sample function Of(t) in our con- 
ceptual ensemble {a^(t)} gives rise to a different function 
Cy(t) as indicated by Eq. (2.4).' Furthermore, fluctuations 
in Cy(t) relative to the mean value of Cy are comparable to 
or less than fluctuations in a|*(t) relative to its mean. 
Consequently, we need only consider variations in the right- 
hand side of Eq. (2.5) that are caused by "small" variations 

in Gy relative to the mean Cy of Cy: 

= E{cte . (2.7) 

y y ^ I ^ 


Such variations in p(y|Gf.) may be efficiently represented by 
a Taylor’s series expansion of the right-hand side of 

Eq. (2.5) in the variable about its mean value - l.e., 

y y ’ 


p(y I a^) 


k =0 k! 


(o^-a^ ) 

y y 


k 


( 2 . 8 ) 


where we have used the definition 

^k 


(k) / I N A 
P ^yl2f^ " 


d(a") 

y 


k 


d 


k 


2 2 
a =a • 

y y 

y^ 


d(o^)^ /2 tt 


2 a| 


y y 


(2.9a) 


(2.9b) 
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according to Eq. (2.5) a and where 




a^=a^ 

y y 


( 2 . 10 ) 


Let us denote the central moments of by 


y 


y y y 


(2.11a) 


= j [E{y^ (t ) I a^} - 


p(Cf) da^ 


(2.11b) 


where, in going to the second line, we have used Eq. (2.4) 
and the vector notation for af(u). If we now substitute 
Eq. (2.8) into Eq. (2.6a), then interchange orders of inte- 
gration and summation, and compare the resulting expression 
with Eq. (2.11), we see that our series expansion of p(y) 
may be expressed as 


P(y) = I 

k=0 


^ 

k! 


p (y I ^f ) 


( 2 . 12 ) 


where P^^^(y|uf) is defined by Eq, (2.9) and is the kth 


central moment of a-t as defined by Eq . (2.11). 


y 


Equation (2.12) is the desired series expansion of p(y) 
Since = Ij the first term in the right-hand side of 


Eq. (2.12) is p(y | q „) 


2 2 
a -a 


, whereas the term corresponding 


y y(i) 

to k = 1 is zero because u ^2 


= 0. Consequently, the first 


two nonvanishing terms of the right-hand side of Eq. (2.12) 
are 


33 



p(y) - p(y |a^) 


+ 


(2.13) 


2 2 
a =a 

y y 


2 


1 ^(2) 

- Pj,2 

y 




CTf ) 


Our motivation for taking Oy = Oy as the expansion point 
of our Taylor’s series is threefold: in this case (1) when 

Cy is a constant 3 the first term is in exact agreement with 
the known Gaussian result, (11) the term k = 1 in Eq. (2.12) 

vanishes identically since =0, and (ill) the second 

^y — 

moment of Cy is minimized about the expansion point = a^; 
hence, the first correction term to the term k = 0 Is^ ^ 
minimized . 


DisQUss-ion. Prom Eq, ( 2 . 5)3 we see that the first term 
in the right-hand side of Eq. (2.13) is the Gaussian density 

function with variance Cy = Cy . Thus, in those situations 
where a-y is a constant — which occur when Of is a constant — 
Eq. ( 2 , 13 ) reduces to the known Gaussian result for stationary 
Gaussian aircraft excitations. In Appendix A, it is shown 
that the low-order correction terms to the Gaussian first 
term in Eq. (2.12) are 


p*^^^(y 


2f ^ 


p(y |o^) 




1 



a 

y 


(2.1^a) 


p(y ) 
2q2 

y 




(2.l4b) 
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(2.15a) 


( 2 ') . P(y|cr.) 

Of.) = 

y 


y 


6^+3 

a" 

y 


a^=a^ 

y y 


p(y |a^) 

Hop^ 


H. 


y 


=v 


y 1 y 


(2.15b) 


and 




P (y I op 

8(a")' 

y 


(a^)' 

y 


- 15 


(a")" 

y 


+ 45 y 15 


y 


2 2 
a =a 

y y 


p(y |a^) 

8(a")' 

y 


H 


X 
6\ a 


. =V^ 

y y 


(2.l6a) 

(2.16b) 


The general form of the above three terms is 


(k) 


(y|c7.) = 


p(y I o^) 


H. 


(2^pk 2k ya 


y; 


V 


y= y 


(2.17) 


where, in Eqs. (2.l4) through (2.17) ^ denotes the 

Hermite polynomial of degree 2k as defined on p. 133 of Ref. 8. 

It may be shown that (at least through the term k = 3) 
the series expansion of Eq, (2,12) is identical to the Gram- 
Charller series of type A — e.g.. Ref. 8, pp . 222, 223. 

However, our derivation of Eq. (2.12) is entirely different 
from the usual derivation of the Gram-Charlier series. Our 
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derivation was based on a Taylov^s series expansion of the 
Gaussian density of Eq. (2.5) where derivates were taken with 
respect to the variance of Cy at a generic point y. Con- 
sequently^ the expansion functions of Eq. (2.9) have the 
property that for sufficiently small low-order central moments 
(k) 

j k ~ 2 ^ 3 ^ ... y a small number of terms in the series of 

Eq. (2.12) should provide an excellent approximation to p(y) 
for every value of y in the range -co<y<co. That is ^ the forms 

(k) 

of the expansion functions p (y\oj^)y k = 2^3^..., (considered 
as functions of y) are the optimal choices for a good fit over 
the entire range -oo<y<oo of interest of the variable y. A pro- 
perty of this type cannot be inferred from common derivations 
of the Gram-Charlier — e.g.^ Ref. 8^ pp . 222^ 223 — and,, in 

factj may be valid only for linear transformations of the class 
of intensity modulated Gaussian processes being considered 
herein — of., von Mises {,32, p. 137~\. 


Example of iwo-ierm Expansion of Response 
Probability Density Functions 

In order to evaluate the terms through, say, k = N in the 

expansion of Eq. (2.8), we require Cy = E{a|} and the central 
moments defined by Eq. (2.11a) through k = N. The quantity 
E{Oy} can be evaluated by integrating over -co<f<oo the power 
spectral density of the aircraft response. In the next 
section of this report, a new method will be described for 

( 2 ) 

evaluating the coefficient y 2 of the first and most 

Important correction term in ?he expansion. The accuracy of 
the two terms of the expansion described by Eq. (2.13) there- 
fore is of considerable interest. 

To ascertain typical accuracy that can be expected from 
the two-term expansion described by Eq. (2.13), we consider 
a one-dimensional analog of Eq. (2.6a), which, we recall is 
the exact expression for p(y) that Eq. (2.13) approximates. 

To relate this one-dimensional analog to the infinite dimen- 
sional integration in Eq. (2.6a), we recall that when fluctua- 
tions in the process aj.(t) are negligible over time intervals 
comparable with the nominal duration of the aircraft impulse 
response h(t), the response process is locally stationary 
and Gaussian {19'] with instantaneous time-varying variance 

Oy E Oy(t) ~ E{y^ (t) |a^(t) } . (2.l8) 
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Equation (2.l8) differs from Eq. (2.4) in that in the quasi- 
statlonary approximation involved in Eq. (2.l8), at each value 
of tjCyCt) depends on the value of the function af(t) ob- 
served only at the same Instant of time — rather than on 
CfCu) observed over -«»<u<^t as in Eq. (2.4). Thus^ Eq. (2.6a) 
may be replaced in this quasl-stat lonary approximation by 


p(y) = 


p(y|o^) p(a^) da^ 


_ XI. 

Q 2 

-- e p(c7^) da. 


0 /2Tra^ 

y 


/ 2^^2 


y 

2ay / 2 ^ 2 

e •^p(a^)da^ 

y y 


:TTa 


y 


(2.19a) 


(2.19b) 


(2.19c) 


where,* in going to Eq. (2.19b) we have substituted Eq. (2.5) 
where Oy is a function of ap as indicated by Eq . (2.l8), and 

in going to Eq. (2.19c) the probability density p(cf) has 
been mapped into the probability density p(ay) with the trans- 
formation between and Cy described by Eq. (2.18). In this 

regard, we note that Eq . (2.l8) implies the existence of a 
generally different value of af for each different value of 
Of — l.e., defines ay as a function of ap. Therefore, the 
probability density p(ap) Implicitly defines from the function 
Cy(af) the probability density p(Cy) shown in Eq. (2.19c). 


V/hen the aircraft excitation is a stationary Gaussian 
process, af is a constant; therefore, from Eq. (2.18) Cy also 
is a constant and p(Cy) is a delta function located at the 
correct value of ay. For this limiting case, the Integration 
in Eq. (2.19c) will yield a Gaussian proabllity density with 
the correct variance. 


*In Eqs. (2.19a) to (2.19c) and the discussion following them, 
we denote the probability density functions of the random 
variables ap and Oy by p(ap) and p(ay), respectively. There- 
fore, p(ap) and p(cy) are different functions of their 
arguments. To keep'^the notational problem from getting out 
of control, we shall generally follow this practice in the 
following pages. That is, probability density functions of 
different random variables will be denoted by p(*) with the 
arguments being the random variables the densities describe. 
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In order to ascertain the typical accuracy that can be 
expected from the two-term expansion of Eq . ( 2 . 13)3 con- 
sider for the probability density p(c^) in Eq. (2.19c) the 
gamma density — e.g., pp, 220-221 of Ref. 33: 


o^r(y) 

0 



ay > 0 


< 0 . 

y 


( 2 . 20 ) 


The gamma density function generally is expressed as a function 
of two "free parameters" which are uniquely determined by the 
mean and variance of the distribution. Instead, we have 

written Eq. (2.20) directly in terms of the mean Oy of the 
distribution and one additional free parameter. This re- 
maining free parameter is the reciprocal of the relative 
variance of Cy — l.e., 


Y 



E{ (a^-a^ ) 

y y 


} 


(2.21a) 


1 


Relative variance of a^ 

y 


(2.21b) 


Thus, when y ^ , the density described by Eq. (2.20) 

approaches a delta function located at Oy = a^. On the other 
hand, when y ~ Is Eq. (2.20) describes the exponential proba- 
bility density. For large finite values of y, it can be shown 
that Eq. (2.20) approaches a Gaussian density in the neighbor- 
hood of . 

y y 
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In Fig. 7a the gamma probability density is shown plotted 
in the usual manner for several values of the parameter y — 
see, for example, pp . 248 and 404 of Ref. 32 . However, we 
have used our notation of Eq. (2.20) in describing the abscissa 
and ordinate in Pig. 7. When plotted as a function of the 

normalized variable of Cy/ay that is "natural" to the gamma 
density as described by Eq. (2.20), we obtain the behavior 
shown in Pig. 8. Notice that each density function shown in 
Fig. 8 has the same mean value, and that Pig. 8 shows the 
relative variance 1/y shrinking as y Increases — as we have 
described above. In the limit y oo p(a|-) approaches a Dirac 

delta function located at Cy = Cy. Figure 8 displays Eq. 

(2.20) in the form relevant to the present work — the plots 
shown in Fig. 7 are included only for comparison with the 
gamma density as it is usually shown. From Fig. 8, we can see 
that the gamma density function of Eq. (2.20) encompasses a 
nice range of shapes to model the probability density of Oy 
for purposes of studying the accuracy of the two-term expan- 
sion of p(y) given by Eq. (2.13). Fortunately, the exact 
density function p(y) also can be evaluated in closed form 
when the gamma density of Eq. (2.20) is substituted into 
Eq. (2.19c) and the integration carried out. 

When Eqs. (2.19c) and (2.20) are combined, it is shown 
in Appendix B that we may express the resulting probability 
density in terms of the normalized response variable 

n ^ (2.22) 



as 

/2y/ TT 

= (/27|n|)^ K ;,(/2 t|ti1) , (2.23) 

2^ ^r(Y) ^ ^ 
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where Kj^C*) is the modified Bessel function of the second kind 
or order n. For n = 3/2^ 5/23... , these Bessel functions 

can be expressed in closed form in terms of elementary func- 
tions. Therefore, whenever y is a positive integer, it follows 
from Eq. (2.23) that py(ri) can be expressed in closed form in 
terms of elementary functions. for n = 1/2, 3/2, 

5/2,..., 19/2 are developed in Appendix C along with expres- 
sions for p^(ri) for y = 1,2, 4, 8, and 9* 

Plots of the probability density p(y/ Cy) given by 
Eq. (2.23) and evaluated from the results in Appendix C are 
shown in Figs. 9 to 12 for values of y of 1,2,4, and 8 
respectively. To compare these results with those produced 
by the two-term expansion of Eq. (2.13), we note first from 
Eqs. (2.5), (2.13), and (2.15), that our two-term approxima- 
tion to p(y) can be expressed as 


2 


P(y) 


2^y 


V 


u(2) 

a + ^ 


2710 ^ 


8(ap^ 


(a")" 

y 



(2.24) 


where from Eqs. (2.11a) and (2.21a), we see that the coeffi- 
cient to the correction term may be expressed in terms of y: 


y 




1 

Y 


(2.25) 


Hence, when we Introduce the normalized response variable of 
Eq. (2.22) our two-term series approximation of Eqs. (2.13) 
and (2.24) becomes 


p(n ) = 





1 + 


1_ 

Hy 


n ‘^-6p ^+3 


( 2 . 26 ) 
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p(i7) =p(y/y^ ) 



FIG. 9. COMPARISON OF EXACT AND APPROXIMATE PROBABILITY 
DENSITY FUNCTIONS, EQS. (2.23) AND (2.26) RE- 
SPECTIVELY, FOR y = 1. DEFINITION OF y IS GIVEN 
IN FIG. 8. 
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p{i7) = p(y/y^ 



v = y//^ 


FIG. 10. COMPARISON OF EXACT AND APPROXIMATE PROBABILITY 
DENSITY FUNCTIONS, EQS. (2.23) AND (2.26) RE- 
SPECTIVELY, FOR Y = 2. DEFINITION OF y IS GIVEN 
IN FIG. 8. 
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p(i7) = p(y/v/^) 


EXACT DENSITY, 
EQ.(2.23) 

APPROXIMATE 
DENSITY, EQ.(2.26) 

GAUSSIAN DENSITY 




- 3.0 - 2.0 


3.0 4.0 


V = )f//^ 


FIG. 11. COMPARISON OF EXACT AND APPROXIMATE PROBABILITY 
DENSITY FUNCTIONS, EQS. (2.23) AND (2.26) RE- 
SPECTIVELY, FOR Y = 4. DEFINITION OF y IS GIVEN 
IN FIG. 8. 


45 




p(i7) = p(y/y^) 


EXACT AND APPROX. 
DENSITIES, EQS.(2.23) 
AND (2.26) 

GAUSSIAN DENSITY 



- 3.0 - 2.0 




FIG. 12. COMPARISON OF EXACT AND APPROXIMATE PROBABILITY 
DENSITY FUNCTIONS, EQS. (2.23) AND (2.26) RE- 
SPECTIVELY, FOR Y = 8. DEFINITION OF y IS GIVEN 
IN FIG. 8. 
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Plots of Eq. ( 2 . 26 ) are shown in Pigs. 9 to 12 for comparl 
son with the exact result given by Eq. (2.23). Also shown in 
each figure is the Gaussian approximation to the density 
function — l.e .3 the first term in Eq. (2.26). Examination 
of Figs. 9 to 12 shows that the two-term approximation of 
Eq. ( 2 . 26 ) provides marginal results for y = I 5 good results 
for Y = 2, and for y = 4 and larger, the approximation is 
excellent . 


We can express gamma in 
deV'iat'ion of Cy — also called 
[ 8 , p. 357] — i.e., from Eq. 


terms of the relative standard 
the GO ef fioient of variation 
(2.21) we have 


y 


( 2 . 27 a) 


= coefficient of variation of . (2.27b) 

y 

Consequently 3 we see from Figs, 9 to 12 that when the coeffi- 
cient of variation of the instantaneous mean-square response 
Oy(t) is unity ^ our two-term series approximation Eq. (2.24) 
provides only marginal accuracy ; when the coefficient variation 
is (1//2) = 0. 707 the accuracy is good^ and when the coeffi- 
cient_ variation of the instantaneous mean-square response is 
(l/}/4) =0.5 or less, the error involved in the two-term 
approximation Eq. (2.24) is, for practical purposes, negligible 

We remind the reader that the probability density function 

of o'y/o'y used to generate the results of Figs. 9 to 12 are 
displayed in Fig. 8. Thus, for all but extremely strong 
variations in Cy(t), the two-term expansion, Eq. (2.24) will 
provide excellent results. To illustrate the excellent fit 
of Eq.^ (_2.24) in the tails of the probability density of 

p(y/|/'ay ) the same curves shown in Pigs. 9 to 12 are plotted 
in Pigs. 13 to 16 on seml-logarlthmlc coordinates. 


Relationship for Expansion Coefficient of 
Correction Term to the Gaussian Density 


Let 

for 

required 


us turn now to obtaining a general relationship for 
cases where no locally stationary assumptions are 
other than that described by Eq. (1.8) — which 
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p(7;)e p(y/. 


— EXACT DENSITY. 
EQ.(2.23) 

— APPROXIMATE 
DENSITY, EQ.(2.26) 

••• GAUSSIAN DENSITY 


IG. 13 . COMPARISON OF EXACT AND APPROXIMATE PROBABILITY 
DENSITY FUNCTIONS, EQS. ( 2 . 23 ) AND ( 2 . 26 ) RE- 
SPECTIVELY, FOR Y = 1- DEFINITION OF y IS GIVEN 
IN FIG. 8. 




/ A) d E (ii')d 





A) d = (tc)d 





p(7/) = p(y/ 


EXACT DENSITY. 

EQ.(2.23) 

APPROXIMATE 

DENSITY, EQ.12.26) 

GAUSSIAN DENSITY 



IG. 16. COMPARISON OF EXACT AND APPROXIMATE PROBABILITY 
DENSITY FUNCTIONS, EQS. (2,23) AND (2.26) RE- 
SPECTIVELY, FOR Y = 8. DEFINITION OF y IS GIVEN 
IN FIG. 8. 



depends only on properties of the turbulence. Thus, the 
results to follow will be valid for aircraft impulse response 
functions of arbitrary long duration — and, therefore, for 
supersonic as well as subsonic aircraft. First, we shall 
obtain a general expression for the mean-square aircraft 
re sponse . 

Geneva! expression for mean- square aircraft response . 
Equation (1.11) is the general input-response relationship 
for spectra, where is the two-sided power spectral 

density of the turbulence, H(f) is the aircraft complex fre- 
quency response function, and ’^’y(f) Is the aircraft response 
power spectral density. The aircraft mean-square response 
is obtained by integrating Eq. (1.11) over -oo<f<co; 

— 

E{y2) E = 0 (f) |H(f)P df . (2.28) 

>■ c/ y j W ' ' 


The bar over Cy is required since E{y^} = Cy is the uncondi- 
tioned expected value of Eq. (2.4). 

From the mutual statistical Independence of the processes 
{Wg(t)}, {ap(t)}, and {z(t)} in our model of Eq. (1.2), it 
follows that the processes {wg(t)} and {wj^(t)} are Independent. 
Therefore, the power spectrum of {w(t)} is the sum of the 
power spectra of {Wg(t)} and {w^(t)} — l.e., 

?> (f) = $ (f) + 0 (f) . (2.29) 

w w w „ 

s f 

Substitution of Eq , (2.29) into Eq. (2.28) shows that the 

mean-square aircraft response is the sum of the mean-square 
responses from the slow and fast turbulence components — l.e., 

a^=a^+a^ (2.30) 

y y^ ’ 


$ (f) |H(f)|" df , (2.31) 

S 


E{yn 

where 
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and 


a 


2 



E{a^ (t)}. 


( 2 . 32 ) 


where 



is the conditional mean-square aircraft response 


to the fast turbulence component conditioned on the random 
function Cf(t) as in Eq. (2.4) — l.e.. 


a 


2 

y 


f 


Oy (t) = E{y^(t ) I a^(u) } , 


_oo<u<t 


(2.33) 


Generat expres s'ion for second centrat moment of oondt- 
ttonaZ instantaneous mean-square response , Since the slow 
turbulence component Wg(t) is Independent of aj.(t), it follows 
from Eq. (2.4) that the conditional mean-square response Oy 
may be expressed as 


a 


2 

y 



+ a 


y. 


(2.34) 


where is the uncondit ioned mean-square response to the 

slow component Wg(t) given by Eq. (2.31), and is the con- 

ditional mean-square response to the fast component defined by 
Eq. (2.33). The second central moment of Cy — defined by 
Eq. (2.11) for k = 2 — therefore can be expressed as 


y ^ 2 ^ - E{ (o^-a ^ ) ^ } 
y y 


= E{ [ (a^ +0^ ) - 


y 


y 


f 


(a^ 

y< 


+a" )]"} 

y^ 


(2.35a) 

(2.35b) 


= E{(a^ 


)U 


(2.35c) 


where we have introduced Eq. (2.34) in the second line and 

have used the fact that a?, is not a random variable. We see 

-y s 

( 2 ) 

from Eq . (2.35c) that y 2 depends only on the aircraft 

^y 

response to the fast turbulence component. 
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Furthermore 5 we see from Eq. (2.24) that the overall co- 
efficient of the correction term to the Gaussian approximation 
of the response density function is proportional to 


( 2 ) 

2 






y 


y 


f 


y 


y< 


+a 


y 


f 


(2.36) 


( 2 ) — 

according to Eqs. (2.30) and (2.35c). Hence, y ^ de- 

2 ^y ^ 

creases as the mean-square response at to the slow 

s 

turbulence component is increased. This behavior is in 
agreement with intuition, since the response to the slow 
component Wg(t) is assumed to be Gaussian; therefore, the 
magnitude of the correction to the Gaussian response term 
decreases as the fraction of the response in the Gaussian 
component Increases. 


Series expansion of conditional instantaneous mean- square 
response, A series of developments in Refs. 34, l8, and 19 
has lead to a series expansion for the conditional instantane- 
ous power spectrum of the "fast" component Wf(t) = Cf(t) z(t) 
of our turbulence model. This expansion is given on p. 23 
of Ref. 19 as 


(f,t 


w 


a^) 


N 

= I 

n= 0 


a ( t ) 
n 

n ! 


(f ) 

z 




(2.37) 


f n ) 

where (f) is defined as the nth derivative of the power 

spectral density {z(t)> — i.e.. 


^(n)(f) A ^ 

z ^^n z 


— and where 




(2.38) 


(2.39) 
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The coefficients a^Ct) in Eq. (2.37) may be expressed in terms 
of the derivatives of Qf(t) by 


= 


1 


n 


(-i47r)^ k=0 


n 

I 


(- 1 ) 


k/n 


d^Q^(t) d^ ^a^(t) 


dt 


k 


dt 


n-k 


(2.40) 


where 


n\ _ n ! 

V (n-k)!k! 


(2.41) 


are the binomial coefficients, and where we remind the reader 
that the conditioning on the random intensity Cf(t) Indicated 
by the left-hand side of Eq. (2.37) implies that c^(t) in 
Eq . (2.40) and the aj^(t) are to be treated as known functions. 
From Eq. (2.40), one may show that for odd Integer values of 
n, we have 

an(t) = 0 , n = odd . (2.42) 

Expressions for the remainder term Rjyj+i(f,t) in Eq. (2.37) are 
given by Eqs. (4.7) and (4.13) on pp. 27 and 29 of Ref. 18. 

The first two nonvanishing terms aj^(t) can be expressed as 

a,(t) = a^(t) (2.43) 

and 

d^ina^(t) 

a,(t) = - aUt) i (2.44a) 

" 8tt^ ^ dt^ 


d^£n[aj(t)] 

^ aj(t) i . (2.44b) 

16 tt^ dt^ 

Hence, we may express the first two nonvanishing terms of the 
series in Eq. (2.37) as 
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w. 




oiit) 


$ (f) - 
z 


327T 


£na^(t ) 
dt^ 




+... . ( 2 . 45 ) 


According to p. 26 of Ref. 19, we can express the condi- 
tional Instantaneous spectrum of the aircraft response to the 
"fast" turbulence component w^(t) as 



( f , t I a^) 


j 0^^(f ,t-T |a^) «>^(f,T) dT , 


(2.46) 


where <|)]^(f,t) is the Instantaneous power spectrum \_S4~\ of the 
aircraft unit-impulse response of Interest : 


$j,(f 


t ) 


h( t--^) 


h( t+J) 


-12'irf T 
e 


dT 


(2.47) 


Substitution of Eq. (2.46) into Eq. (2.37) gives — Ignoring 
the remainder term ~ 


^ (f ,t 
^f 




N 

= I 

n= 0 


1 

n! 


a ( t-T ) 
n 


z 


^j,(f,T) 


dx 


(2.48) 


To obtain an expression for 
mean-square response to the 
Eq. (2.33), we Integrate <|) 


the conditional Instantaneous 
fast component Wf(t) as in 
(f,t|a^) over all f: 
f 


E (t) = E{y^(t ) I a^(u) } 

^ ^ f ^ 


-oo<u< t 




df 


or 
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at = (yt (t) = I ^ j a^(t-T)| $j^(f,T)df dT 


y. 


n=0 

N 1 

= ^ J 

n=0 J 


n 


J 

_oo 


a (t-x) Yv, (t^) dT , 
n 'hyZ ,n ’ 


(2,^9) 


where, in going to the third line, we have used Eq. (2.48) and 
have Interchanged the order of Integration, and in going to 
the last line, we have used the definition 


Yk ( t ) = 

' h, z , n 


(f ) 

z 


df 


(2.50) 


which is a function of the nth derivative of the power spectral 
density <l>z(f) of the turbulence component {z(t)> and the in- 
stantaneous spectral density cE>j^(f,t) of the aircraft unit- 
impulse response function h(t). The result of Eq. (2.49) 
Ignores the contribution of the remainder term RI^J 4 .l(f,t) in 
Eq . (2.37). Equation (2.49) is the desired series expansion 
for the conditional instantaneous mean-square aircraft response 
to the "fast" turbulence component Wf(t) = ap(t) z(t). 


Series expansion of second moment of conditional instan- 
taneous mean-square response . An expression for the coeffi- 
( 2 ) / — 

dent u ^ -2 / ( o ^ )^ of the correction term for the Gaussian 

aj/ y 

probability density in Eq . (2.24) is given by Eq. (2.36), 
where we see that the second central moment of a^ 


yf 


is required 


We may obtain a series expansion for the second moment of a^ 

y f 

by taking the expected value of the square of Eq. (2.49): 


E{ (a^ 

y 


)M = 


f 


n ^=0 n 2=0 


CO 


n^ In^ ! j 


E{a (t’-T )a (t’-T^)} 
n i n 2 

1 2 


X Y (t ) Yi„ (t ) dT , dT 

’^h,z,n^ 2' 1 


(2.51) 
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where the expectation is taken with respect to fluctuations 
in the random function ap(t), and where t has been replaced by 
t’. ¥e are primarily interested in situations where the 
random process {Cj^Ct)} is stationary. In these cases, we may 
take t = t^ " hence. 


E{a^ (t'-T^) a^^ (t’-T^)} 


E{a (t ) a (t+T -T^ )} 
n n 1 2 

1 2 


(2.52) 


which is Independent of t . 

Intro duo t-ion of loaatZy stationary assumption , Our main 
interest in these developments is situations where the locally 
stationary assumption of Eq. (1.8) is valid. Prom Eqs. (1.8a) 
and ( 2 . 43 ) through (2.45), we see that the locally stationary 
assumption permits us to include only the term n = 0 in 
Eq. ( 2 . 49 ) -i.e.. 


a 


2 

y 


f 


(t) 


E{y^(t ) I a^(u) } 


_oo<u<t 


a^(t-T) ^(t) dT , 


(2.53) 


where we have used Eq. (2.43) and the notation 


0 00 




5>^(f) <5^(f,t) df 


(2.54) 


which follows from Eq. (2.50). The expected value of Cy^ 
[with respect to fluctuations in a^(t)] is 


a 


2 



E{a^ (t)} 
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hence 
Eqs . 


, substituting the mean value of 
(2.53) and the above relationship 



(t ) , 


we have from 


o 


2 



(t) 



[a^( t-T ) 





dT 


(2.55) 


According to Eqs. (2.24) and (2.35c) we require the 

( 2 ) 

second central moment y'a which is given by the mean-square 

value of Eq. (2.55). However, the right-hand side of Eq. 

(2.55) is the convolution of [af(t) - Of'] and “ i.e., 

Eq. (2.55) may be thought of as describing the response of a 
linear two-terminal time-invariant system with unit-impulse 


response yj^ z*^^^ stochastic input 

we may express the power spectral density 


cess 

slth 


{(j| (t) - Oy } as the product of 
^>^^2ff) of {a|(t) - a|-} and | y^^ ^ 


the 

(f) I 


[af(t) - Of]. Thus, 
^^2 (f) of the pro- 

yf 

power spectral den- 
^ , where 


. . A 




i27TVt 

e 


dt 


(2.56) 


is the (Inverse) Fourier transform of the ’’system impulse 
response” y, „(t): 


$ 



(f) 


$ 2 
‘’f 


(f) 




(2.57) 


The me an- square value of {of^ (t) - o^ } 

Yf yf 

Eq. (2.35c) is the second central moment 
integrating Eq. (2.57) over all f: 


which 

^y 


according to 
is obtained by 


n^2) 

y 


^>^2(f) 

Of 




df 


(2.58) 
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Equations (2.28) and (2.68) are the required retation- 

( 2 ) 

ships for evaluation of the coefficient y 2 of the 

the correction term that occurs in Eq. (2.24) to the Gaussian 
probability density function of the aircraft response y(t). 
The turbulence characterizations required for evaluation of 

y 


and n ^2 cere the power spectral density ^^^(f) of the tur- 


y 


bulence velocity [w(t)} and the power spectral density ^^ 2 (f) 

of lo^( t) - where Oj'( t) is the intensity modulation ^ 

in the fast component w^(t) given by Eq. (1.2). Methods for 
computing ^ ^^ ( f ) from turbulence recordings are described on 

pp . 79-82 of^Ref. 19. The quantity defined by 

Eqs. (2.54) and (2.56). Notice that y^ ^ ^ depends both on 

the power spectral density ^^(f) of the turbulence component 
{z(t)] and on the aircraft unit-impulse response function as 
may be seen from Eq. (2.47). Thus ^ y, (f) or y, (t) char- 

fZ j ^ /2 j ^ 

acterizes the aircraft impulse response with respect to the 
turbulence component {z(t)}. The locally stationary assump- 
tion of Eq. (1.8) has been used in deriving Eq. (2.58). 


Alternative form for system characterization. Generaliz- 
ing the definition of Eq. ( 2 . 56 ) to include y^ _ — which 
is defined by Eq. (2.50) — we have 


h , z ,n 


/ ^ A 

^h, z 5 n ^ 


/ , ^ i27TVt ,, 


I I df 


rh">(f) 

J Z 

.00 

J e^2wvt 

-^00 

_oo 


i27TVt 
e dt 


df 


(2.59) 


However, using Eq. (35b) on p. 29 of Ref. 3^: 
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oo 






12Trvt 


dt 


= H(f-j) H*(f+|-) 


(2.60) 


where H(f) Is the complex frequency response defined by 
Eq . ( 1 . 9)5 have from Eqs . (2.59) and (2.60) 


^h , z 5 n ^ ^ 


■ / \ 

$^^^(f) H(f-^) H*(f+^) df. 

^2 c. cL 


( 2 . 61 ) 


The special case of Eq. ( 2 . 6 l) for use in Eq. (2.58) is 


y 


h, z 


(v) 


r ^ 

H(f-j) H*(f+|-) df , 

_ OO 


( 2 . 62 ) 


which is the desired alternative form of y, (v) for use in 

' h , z 

Eq . ( 2 . 58 ). Equation (2.62) may be Interpreted with the aid 

of the material contained in Sec. 4 of Ref. 34. 


Limiting Cases of Expansion Coefficients of the 
NonGaussian Term 

Limiting cases of the central result, Eq. (2.58), provide 
insight into that result. Three limiting cases are discussed 
below . 

Case 1: w^(t) vs stationary and Gaussian. Since the 

process {z(t)} in our model of Eq. (1.3) is stationary and 
Gaussian, the "fast” turbulence component {wf(t)> is stationary 
and Gaussian when 

c^£.(t) = constant = . (2.63) 

In this case, the power spectrum ^ 2 (f) of the random process 

2 — Y ^ f 

{of(t) - Of} is zero. Therefore, according to Eq. (2.58), 
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( 2 .) 

is zero, and according to Eq . (2.24), the probability 
density of the aircraft response y(t) is Gaussian. 

Case 2: power spectrum of {z(t)} is white* This case 
is of little direct Interest for the representat ion of tur- 
bulence since {z(t)} would normally be taken to have the 
appropriate (transverse or longitudinal) von Karman spectral 
form. Nevertheless, it is of general Interest as a limiting 
case in the study of Eq. (2.58). For this case, we have 


(f) = constant = $ (•) 


Using Eq. (2.64), we have from Eq. (2.54) 


(2.64) 


= $^(*) h^(t) 


(2.65) 


where, in going to the second line, we have used Eq. (12a) on 
p. 26 of Ref. 3^ applied to the aircraft unit-impulse response 
h(t). Therefore, according to Eq. (2.56), we have for the 
present case 






hMt) 


-.00 


12'ITVt 

e 


dt 


( 2 . 66 ) 


which is the (Inverse) Fourier transform of h^(t) multiplied 
by the constant spectral density $ ( • ) . 

Equation (2.66) may be used directly in Eq. (2.58) to 
evaluate y 2 . Alternatively, we may define the "autocorrela- 

tlon function" of the square of the (deterministic) unit- 
impulse response h(t) as 




h^ ( t ) h^ (t+T ) dt 


(2.67) 


which is necessarily an even function of t 
theorem — e.g., Eq. (135) on p. 5^ of Ref. 
from Eqs. (2.66) and (2.67), 


Using Weiner ’ s 
34 — we then have 
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it>jj2(T) dx , 


( 2 . 68 ) 


r 2 ) 

which may be used to evaluate y 2 by Eq. (2.58). Notice that 

h^(t) is a nonnegative function. Therefore, for reasonably 
well behaved h(t), the main contribution to yj^ z(v) will be 
in the region about the origin v = 0, and the nominal bandwidth 
of|Yh z(^)|^ will be of the order of the reciprocal of the 
nominal duration of h(t). Therefore, this same low-frequency 
region of the spectrum $^ 2 (v) is relevant in the computation 

of by Eq. (2.58). 

Case 3: fluctuati-ons -in O-r(t) are negZ'igtb'le over the 

duration of h(t). When fluctuations in the intensity modula- 
tion ap(t) occur sufficiently slowly, the power spectral 

density of {a|>(t) - Of} has all of its area concentrated near 
zero frequency. The limiting case in this situation is 


^>^z(v) = E{(a^-ap^} 6(v) , 


(2.69) 


where 6(v) is the Dirac delta function located at v = 0 . Sub- 
stitution of Eq. ( 2 . 69 ) into Eq. ( 2 . 58 ) yields 


E{(o|-ai)n 

y 


5 (v) I Yh z I ^ 


= E{ (a|.-ap^} I ^ 


(2.70) 


However, from Eq. (2.62) it follows directly that 


= I |H(f)y df , 


(2.71) 


63 



which is real; hence, Eqs. (2,70) and (2.71) yield 


y^2 E{ (a^-ap ^ } 

y 


|H(f) 



(2.72) 


which is the desired result. 

From Eq, (2.71), we recognize ^(0) as the mean-square 
response of the aircraft to the turhu%ence component {z(t)]. 
When Eq. (2.72) is combined with Eq. (2.24), we see that the 
resulting expression for the probability density function of 
the aircraft response is identical with that given by Eqs. 
(4.50) and (4.51) on p. 48 of Ref. 19. These results were 
derived under the assumptions that (i) fluctuations in Of(t) 
occur slowly in comparison with those of z(t), and (ii) varia- 
tions in oj^(t) are negligible over durations comparable with 
the aircraft impulse response duration. The more general 
result of Eq. (2.58) relaxes, completely, the requirements for 
assumption (ii) . However, assumption (i) — which is described 
by Eq. (1.8) of this report — has been used in obtaining 
Eq. (2.58). We again emphasize that the assumption of Eq. 

(1.8) depends on turbulence properties alone and is believed 
to be generally satisfied by atmospheric turbulence . 


Series Representation of Expansion Coefficient 
of the NonGaussian Term 

The limiting case of Eqs. (2.70) — (2.72) suggests a 

( 2 ) 

series expansion for y 2 that may be useful in situations 

^y 

where assumption (11) above is almost or only marginally 
satisfied. Equation (2.70) depends on the frequency content 
ofly, (v) I ^ only at v - 0 . This fact suggests a series 

representation of y 2 ^ obtained by first expanding |y (v)|^ 

Oy n , ^ 

in a Maclaurin series and then integrating term by term. 

The Maclaurin expansion of | y, (v)|^ may be expressed as 

in 5 z 




00 


I 

n=0 


n 

V 

n! 


d 


n 


dv 


n 




2 




v=0. 


(2.73) 



By Leibniz’s rule for the nth derivative of a product (p. Ill 
of Ref. 35)5 we have 


d I ry / . N I 2 _ d. 


n 


dv 




n ' ' h, z 


dv 


n 




n , , d*^yg^^(v) d"-\'^(v) 


Jo (E) 


dv 


k 


dv 


n-k 


(2.74) 


However, from Eq. (2.62), we have 




dv 


v = 0 


> k 

* (f) ^ 


H(f - ^)H*(f + |) 


df 


v=0 


(12v)*" $^(f)m^‘"^(f)df. 


(2.75) 


where we have used Eq. (79) on p. 269 of Ref. 36 as applied 
to the deterministic (complex) function H(f), and where 

m^^^(f) is the kth power-moment spectrum of the aircraft unit- 
impulse response function h(t). Properties of power-moments 
spectra are discussed on pp . 26^-269 and 281-288 of Ref. 36 . 
Their name arises from the fact that their Integrals over 
_oo<f<oo satisfy 


^00 . -00 

m^^^(f)df = t^h^(t)dt, k=0,l,2,..., (2.76) 


where the right-hand side of Eq. (2.76) may be Interpreted as 
the kth time-moment of the Instantaneous "power" h^(t) of 

h(t). Furthermore, the power-moments spectra m^^^(f) are real 
and even functions of the frequency variable f. Prom Eq. 
(2.75), we also have 
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dv 


k 


v=0 


(- 1277 )^ 


$ (f (f )df . 
z h 


(2.77) 


Let us define 


^00 


( 2 . 78 ) 


Then combining Eqs. (2.74), (2.76), (2.77), and (2.78), It 
follows that 


_d 

dv 




v=0 


(12.)" I (-l)k (") 

k=0 \ / 


^(k) ^(n-k) 
h 3 z h 3 z 


(2.79) 


Hence, from Eqs. (2.58), (2.73), and (2.79), we have 


f O \ oo ^oo 

= I ^ (12trv)"$ ,(v)dv 

y n=0 • f 


1 

^ I (-1)'^ 

k=0 



„(k) (n-k) 

h, z n , z 


( 2 . 80 ) 


However, the aut ocorrelat ion function <p 2 (t) of the random 

^ ^f 

process {ap(t) - Qf>} is the Inverse Fourier transform of the 
power spectrum $ 2 (v) — i.e., 

Of 

f oo 

(v)e^^^^^ dv . (2.8l) 

f 


Differentiating both sides of Eq. (2.8l) n times with respect 
to t yields 
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4) 2(t) 
dt 


(i2TTv)^ 4> 2(v)e^^^^^dv; 


( 2 . 82 ) 


hence , 


<t>o?^(0) 

f 


(12ttv)^ <J>^2('^)dv 


(2.83) 


Let us define 


(n) A 


n 


h 3 z 


= I (-1) 


k /n\ „(k) T.(n-k) 


k=0 


k/ h, z h, z 


(2.84) 


Combining Eqs. (2.80), (2.83), and (2.84) we have, finally. 


/f) = ^ 4k?ko) 

y n =0 f 


’h , z 


(2.85) 


( 2 ) 

which Is the desired series expansion of y 2 • 

The power-moments spectra are defined on p. 260 of Ref. 

36 as 


m^ C f ) 


t^$^(f,t)dt, k= 0 ,l, 2 . 


( 2 . 86 ) 


Hence, by forming the kth moment of ^(t) and using Eq . 

( k ) -i-i j 

( 2 . 5 ^), we see that ^ and y, (t) are related by 

h , z n , z 




(f) [ t^$^ (f ,t )dtdf = , (2.87) 

h h , z 
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where Eqs . (2.86) and (2.78) have been used In the second 

equality. According to Eq . (2.87), the kth moment 

with respect to time of 

Discussion , The terms in Eq . (2.85) for odd values of n 

( n ) 

are identically zero. To show this, we show that s, =0 

^ h, z 

for n = odd. Consider Eq. (2.84) and let k ' = n-k; hence, 
k = n-k’. Therefore, 


n ! 


n ! 


(n-k)!k! (k’)!(n-k')! 


n 

k’ 


( 2 . 88 ) 


Consequently, we may express the summation in Eq. (2.84) as 


n 


,(n) ^ V 
’h, z ^ 


(- 1 ) 


n-k’ /n \ „ (n-k ’ ) ^ (k ’ ) 


k’=0 


k 


r 


h , z 


h , z ^ 


(- 1 )" I (- 1 ) 

k'=0 


n 


k ’ /n \ p (k ' ) p (n-k ’ )' 
k ’/ h , z h , z 


(2.89) 


Comparison of Eqs. (2.84) and (2.89) shows that 


h, z h, z 


( 2 . 90 ) 


from which it follows that 


4"^ = 0 

h, z 


n=l,3,5,... . 


(2.91) 


Therefore, the odd numbered terms in the summation of Eq . 
(2.85) are Identically zero. 

It is instructive to consider the first two nonvanishing 
terms of Eq. (2.85). For n = 0, it follows directly from the 
definition of the autocorrelation function that 


<)>^F(0) s 
f 


<l>^2(0) = E{[a|(t)-apn 


(2.92) 
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Furthermore 5 one may show [e.g., p. 21 of Ref. 37] that 


0 



( 0 ) 



(2.93) 


which is the negative of the mean-square velocity of a^(t). 
In addition, we have from Eq. (2.84), ^ 


JO) ^ 1 ^( 0 ) 

h , z h , z 


(2.94a) 


/■“ r 

^,(f) 


f , t ) dtdf 


$^(f) |H(f) I ^df 


(2.94b) 


according to Eq. (2.8?) above and Eq. (12a) on p. 26 of Ref. 
34. Furthermore, from Eq . (2.84) we also have 


(2) ^ (0) (2) _ p/p(l)V . p(2) (0) 

h,z h,z h,z \ h,z/ h,z h,z 



p(2) 

n , z 




(2.95) 


According to Eqs. (2.91), (2.94a), and (2.95), the first non- 
zero terms of Eq . (2.85) therefore can be expressed as 


( 2 ) 

U 2 
a 

y 


\ h , z 


+ 


<t>h4o) 

f 


,( 2 ) 
h , z 
.( 0 ) 
h , z 



( 2 . 96 ) 


Prom Eqs. (2.92) and (2.94), we see that the first term in 
Eq. ( 2 . 96 ) is identical to the limiting case of given 
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by Eq. (2.72), which applies to situations where variations 
in af(t) 'are negligible, over durations comparable with the 
aircraft impulse response duration. 

We consider, now, the second term in the right-hand 
side of Eq. (2.96). Using Eq. (2.87), we see that 





where 


t = t. 


A —00 






( 2 . 98 ) 


is the t Ime-centrold of z(t). The equivalence of the two 
right-hand sides of Eq. (2.97) is easily proved by expanding 
(t-t)^ in the second line and using Eq. (2.98). Thus, Eq. 
( 2 . 97 ) represents the second central moment of the normalized 
"mass density" which is analogous to the standard 

deviation of the "density function" y, (t). 

It is shown on pp, 100-101 of Ref. 19 [Eq. (A. 18 ) in 
particular] that the quantity 2 ^ ( 0 )/(f) ^ 2 ^ ( 0 ) is about 

Of 

one-third of the nominal correlation Interval of the process 
{cr^(t) - a^} , Consequently, we have 
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r„(2) 

f h,z 

Cj)^2^(0) 

L h^z 


"nominal duration" of y, (t) 

n 3 z 

"correlation interval" of {a|.(t) - a^} 


(2.99) 


where, here, the sign t should be read as "is of the order of," 
and where we have considered the "nominal duration" of Yh 
to be about three times the duration of Y]^ z^^) measured by 
the square-root of its second central moment. Consequently, 
when the correlation interval of the random process 

{cr^(f) - a^} is large in comparison with the nominal duration 

of Yh z ("*^)3 second term in the right-hand note of Eq. 

’ ( 2 ) 

( 2 . 96 ) is negligible and the approximation to y 2 given by 

*^y 

Eq. ( 2 . 72 ) is adequate. Moreover^ we see from Eqs. (2.72), 
(2.96) j and (2.99)^ that the approx-Cmat-ion Eq. (2.72) wilt 

( 2 ) 

tend to overestimate the positive quantity y 2 . 

^y 

Finally, to obtain a physical understanding of the dura- 
tion of Yh z("*^) consider the case where the power spectrum 
of the process {z(t)} is white — i.e., “ 

constant. In this case, we see from Eq . (2.65) that 




(2.100) 


hence, when ^^(f) ^ constant over the "passband" of h(t), 

the above statements pertaining to the duration of 'Y'h.z(t) 
can be interpreted as pertaining to the duration of h^(t). 


In cases where the mean-square velocity of the process 

2 1 ( 2)1 

{cf(t)} does not exist, 2 ^( 0 )| does not exist as may be 

Qf 

seen from Eq . (2. 93) 5 and the series expansion Eq. (2.85) and 
the two-term approximation Eq. ( 2 . 96 ) are of little use. An 

( 2 ) 

alternative series expansion of y 2 is derived in Appendix D. 
This alternative expansion has the advantage that the first 
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correction term to the limiting case of 


( 2 ) 

given by Eq. 

y 


(2.72) does not require the existence of the mean-square 
velocity of the process {c|*(t)}. The alternative expansion 
derived in Appendix D also is more amenable to calculations 
based on numerical estimation of the power spectrum and auto- 
correlation function of {a^(t)} than the expansion, Eq. (2.85). 


In this section, we have provided a detailed methodology 
for estimating the aircraft response probability density 
function from appropriate characterizations of the turbulence 
excitation. Comparable techniques are developed in Sec. 6 
for threshold mean exceedance rates of the aircraft response. 
For cases where fluctuations in cjp(t) in our turbulence model 
of Eqs. (1.2) and (1.3) are negligible over durations compar- 
able with the duration of the aircraft Impulse response h(t), 
expressions for mean exceedance rates provided on pp . 36 to 
46 of Ref. 19 are applicable. 
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MAXIMUM LIKELIHOOD ESTIMATION OF THE INTEGRAL SCALE AND 
INTENSITY OF von KARMAN TURBULENCE 


Here, we shall derive an optimum method for estimation 
of the Integral scale and Intensity of a well-behaved turbu- 
lence record with no appreciable low-frequency component 
Ws(t) present. The vertical time history shown In Pig. 17 
Illustrates such a record. To develop the method, we shall 
assume that the functional form of the power spectrum of 
the record Is known — e.g., the transverse or longitudinal 
von Karman spectral forms. The problem, then. Is to develop 
an optimal method for estimating the parameters that deter- 
mine the spectrum. In the case of the above two von Karman 
forms, two parameters — the Integral scale L and Intensity 
a — completely determine the spectrum. Several ad hoc 
methods for estimating values of the Integral scale are 
described on pp . 356 to 360 of Ref. 22. 

The method derived below Is based on the Intuitively 
appealing procedure called the method of maximum likelihood ■ 
e.g.. Ref. 38 — which was orglnally Introduced by Gauss. 

This method Is usually treated In the context of optimal 
estimation of the parameters — e.g., mean and variance — of 
a probability density of known functional form from a random 
sample "drawn" from the density. As we show below. It also 
may be used to estimate the parameters L and a of a turbu- 
lence record whose power spectrum Is of known functional 
form. 


Maximum Likelihood Estimation of the Parameter 
in a Probability Density Function 

Let us first briefly review the method In Its usual 
context. Let p(x| 0 ) be the probability density of a random 
variable x, where 0 Is a parameter In the density function, 
such as the mean value. Let X^,X 2 ,...,X be n samples drawn 
from the population whose probability density Is governed 
by p(x|0). The likelihood function Is defined as 


L(x^,X2,...,x^|e) 


n 

n p(x,|e) 

j=i 


(3.1) 
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which Is a function of the n sample values Xj 5 X 2 , . . . j,Xj^ and 
0 , where the functional form p(x| 0 ), considered as a function 
of X and 0, is assumed known. The problem Is to determine 
the "best estimate" of the unknown parameter 0 from the n 
sample values Xj ^X^ , . . . ,X^. When 0 Is assumed known, the 
right-hand side of Eq. (3.1) Is the Joint probability density 
of the sample Xi ,X 2 , . . . ,Xj^. The method of maximum likelihood 
assumes that when 0 is unknown, the "most likely" value of 0 
is the value that maximizes the probability of the observed 
sample. That Is, when Xi,X 2 ,... 5 Xn are substituted Into the 
right-hand side of Eq. ( 3 .D 5 the "most likely" value of 0 
Is the value that yields a maximum value of L(Xi ,X 2 , . . . jX^ | 0 ) . 
Generally, It Is most convenient to maximize log 
L(Xi ,X 2 5 . . . ,Xn I 0 ) rather than L Itself. Thus, the "most 
likely" value of 0 Is obtained by solving the equation 

9L(Xi,X2,. . . ,X^|0) 

= 0 . (3.2) 


When more than one solution to Eq. (3-2) exists, the solution 
chosen Is the one that maximizes L when substituted back 
Into L(X ,X ,...,Xj^ ). Equation (3-2) Is referred to as a 

likelihood equation. For Eq . (3.2) to yield a maximum, it 

Is necessary that (9^L/99^)<0 at the value of 0 determined 
by Eq. (3- 2) . 


Joint Probability Density of Unsmoothed Turbulence Spectra 

Let us turn now to the problem of estimating the in- 
tegral scale and Intensity from a recorded turbulence time 
history such as the vertical record shown in Fig. 17 which 
exhibits no discernible slow component w^(t). To begin 
with, we shall assume that the record Is^a time history 
of finite duration drawn from a stationary Gaussian process; 
later, we shall argue that the final result Is not parti- 
cularly sensitive to the stationary assumption. Furthermore, 
we shall assume that the duration T of the record is suf- 
ficiently large so that no appreciable bias distortion of 
its spectrum Is caused by operating only with a finite seg- 
ment of T secs duration. That is, we shall assume that for 
appropriate (unknown) choices of Integral scale and intensity, 
the expected value of the power spectrum of our finite seg- 
ment of duration T Is equal to the actual power spectrum of 
a record of infinite duration. 
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FIG. 17. LOW-ALTITUDE TURBULENCE RECORDS UNDER CONVECTIVE 

CONDITIONS. [AIRCRAFT SPEED 129 m/sec (422 ft/sec).] 
(Ref. 23, Fig. 4, p. 282.) 
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The problem^ then, is to estimate L and from a segment 
T secs long drawn from a stationary Gaussian process. The 
power spectral density of the process Is of known functional 
form, and is determined by the two parameters L and . 

Our interest, of course, is to determine the values of L and 
for the underlying process from which our finite sample 
has been drawn. Thus, we must consider the statistical 
properties of a conceptual ensemble of stationary Gaussian 
segments T secs long that are drawn from the underlying 
process . 

From each such segment, we first generate, conceptually, 
a periodic function by repeating our T sec segment end-on-end. 
That is, if we denote a typical Gaussian segment by w(t), 

T T 

- then from each such segment we generate a periodic 

function satisfying 

w(t+pT) = w(t) , p=0,±l,±2, . . . . (3.3) 


Since each such function in our new ensemble Is periodic with 
period T, we may consider the statistics of Its complex 
Fourier series coefficients c^^, n=0 , ±1 , ±2 , . . . which are the 
complex amplitudes of the Fourier series components occurring 
at frequencies of f=±n/T, n=0,l,2,... — l.e.. 


w( t ) = y c e 

L 

n=-oo 


12TTnt/T 


where 


n T 


T/2 


-T/2 


/ , ^ -127T nt/T . , 

w( t ) e dt 


= a -lb 
n n 


(3.4) 


(3.5a) 


(3.5b) 


where aj^ and bj^ are real. First, we note that since the 
original process w(t) is Gaussian, the complex Fourier co- 
efficients c^ must be jointly Gaussian complex variables 
since the operation on w(t) described by Eq. (3.5a) Is a 
linear transformation. In particular, for a given value of 
n, the probability density of the real and Imaginary parts 
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and bj^ of the complex coefficients Cn must be governed by 
a Joint normal probability density on the phase plane. 
Furthermore^ since the original process w(t) is stationary , 
there can be no preference between the cosine and sine com- 
ponents aj^ and b^; that is, for every harmonic n, the sta- 
tionary property of w(t) requires that the Joint probability 
density of the cosine and sine coefficients on the phase 
plane be rotationallv invariant. It follows that a^^ and b^^ 
must be uncorrelated and that each is governed by a normal 
probability density with the same mean value of zero and 
variance of, say, o^/2 . For a discussion of the bivariate 
normal density, see, for example, pp . l47 and 1^8 of Daven- 
port and Root \_Z9~\, 

Consider, now, the Fourier transform of a typical tur- 
bulence segment — which we now defined to be zero outside 
the Interval (-T/2 )^t< (T/2 ) . From the sampling theorem in 
the frequency domain — e.g., p. 33 of Woodward [4C>] — we 
see that at the values of f^ = ±n/T described above, the 
Fourier transform of our turbulence segment is T times the 
values of the complex Fourier series coefficients c^^. This 
fact also is Immediately evident from Eq. (3- 5a). Further- 
more, we see from the same sampling theorem that the Fourier 
transform of our truncated segment is completely determined 
by the complex coefficients Cj^, n=0 , ±1 , ±2 , . . . . 

Consider, now, the statistical properties of an estimate 
of the power spectral density of the turbulence process, where 
the estimate is the so-called perlodogram defined as 


S(f ) 


k 1 

T 


rT/2 

w( t ) 

J 

-T/2 


-i2TTf t 
e 


dt 


(3.6) 


— e.g., p. 107 of Davenport and Root [35']. At frequencies 
fn = n/T, S(f) is related to the Fourier series coefficients 
of Eq. (3.5) by 


S(f ) E 
n 


S = 
n 


Tc 


n 


= T c 


n 


= T[a^+b^] 


n n 


n==l,2, . 


(3.7) 
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L 



According to the above comments ^ an and bj^ are Independent, 
and each is governed by a normal probability density with 
zero mean and variance Pn/2. Hence, the random variable 



+ b^ 


n 


(3.8) 


Is governed by the exponential probability density function 




(3.9) 


as is shown, for example, on p. 53 of Lawson and Uhlenbeck 
[47], who, incidentally, incorrectly refer to the density 
of Eq. ( 3 . 9 ) as a Rayleigh density. It follows from Eqs. 
( 3 . 7 ) to ( 3 . 9 ) that values of also are governed by an 
exponential density function, say. 


, S >0 
’ n— 


, S ^<0 , ( 3 . 10 ) 

where = l/(Tp^). Finally, we note that the random vari- 
ables Sj^, n = 0,1,2,... are mutually independent since the 
real and imaginary parts of the complex Fourier series co- 
efficients Cn, from which the are generated, are uncor- 
related Gaussian variates. A proof of the fact that all 
pairs of random variables a^^, a^y^, n^m; b^, bj^^, n?^m; and aj^, 
bni, all n,m, are uncorrelated for n, m^O is provided in 
Appendix E. Also see p. 119 of Goldman [42]. 

Our original time segment w(t), (-T/2 )£t< (T/2 ) from 
which S(f) is generated by Eq. (3.6) can be reconstructed 
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from its complex Fourier series coefficients c^ = aj^-lbj^, 
n = 031 , 2 ,... . However, we have shown above that for all 

m, n^O, all nonldentlcal pairs of an and bj^ are statistically 
Independent. Furthermore, the random variables an and b^^ 
both are governed by a normal probability density with zero 
mean and the same variance (p^/ 2 ) = (l/ 2 XnT), which Is 
uniquely determined by n = 0,1,2,... and T. Consequently, 

all statistical Information about a segment of w(t), 

(-T/2 )£t< (T/2 ) , Is contained In the sequence of coefficients 
Anj n = 0,1,2,... . It follows that the sequence of pro- 
bability densities p(Syi)^ n = given by Eq. (3.20) 

must contain all possible statistical information about the 
periodogram S(f) defined by Eq. (3.6) since it provides a 
complete statistical description of the random process w(t)^ 

( -T/ 2 ) ^t< ( T/ 2 ) j from which S(f) is generated. 

Since all Sj^ are Independent for n>0, the joint pro- 
bability density function of the "vector random variable" 

{ S j , S 2 5 . . . , Sj^} Is 


S^.>0, j = l,2,. . . ,N 


S <0, j=l,2,. . . ,N 

(3.11) 

where we have used Eq . (3-10), and where N can. In principle, 

be Infinity. Finally, we note that Aj Is the reciprocal of 
the expected value of Sj — as Is easily shown In Eq. (3.10): 

A. = i , j=l,2,...,N . ( 3 . 12 ) 

J E{S.} 

J 

Equation (3.11) is the Joint probability density func- 
tion of samples Sj of the periodogram defined by Eq. (3.6) 
at frequency values of 

fj = j/T , J=1,2,...,N , (3.13) 
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where the periodogram samples taken at the frequencies de- 
fined by Eq. (3.13) contain all of the Information in the 
periodogram. Equation (3*13) relates the Xj in Eq . (3.11) 
to the expected values of the unsmoothed power spectral 
density samples Sj = S(fj). 

Likelihood Equations for a General Class 
of Turbulence Spectra 

We may rewrite Eq. (3.11) as 


p(Sj . ,Sjj) = e 


-(AjSj+A 2S2 + . . .+AjjSjj) 


(3.14) 


where, for convenience, we have left out the statement that 
the right-hand side Is valid only for Sj>0, j* = 1,2,...,N. 
The logarithm to the base e of Eq. ( 3 . 1'^T Is 

, . . . ,S ) ] = 


[ £n ( X j X 2 . . . Xj^ ) ] - [ X j S j +X 2 3 ^ + . . . +X^S^] 




[£n(S^S^. . .Sj^)] + 


N^N- 


— + — +. . . + 


'N 


J 


(3.15) 


where. In going to the last line, we have used Eq . (3.12) 

and the notation 


S^. E E{S^.} , J=1,2,...,N . 


(3.16) 


Let us now Introduce a class of power spectrum func- 
tional forms 


E{Sj} E S^. = q^L F^.(L) , j=l,2,..,,N , (3.17) 
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where a Is the rms value of the turbulence, L is the integral 
scale, and for a given frequency index J, Pi (L) is a function 
of L but is independent of a. Equation (3*17) includes the 
(two-sided) von Karman transverse and longitudinal spectra — 
e.g., pp. 83 and 93 of Ref. 18 — 




1+188.75 

[l+70.78L^k^]^^'^® 


(3.18) 


(k) = a^L — 

[l+70.78L^k^]^-^^ 


(3.19) 


where we have used wavenumber k Instead of frequency f as is 
conventional in turbulence work. According to Eqs. (3-17) 
to (3.19)5 for the von Karman transverse spectrum^ we have 


1+188.75 L^k^ 

F (L) = i— 

^ [l+70.78L^k!]“'^® 

J 


( 3 . 20 ) 


and for the von Karman longitudinal spectrum, we have 


P.(L) = £ 

[l+70.78L^k?]^/^ 

tJ 


( 3 . 21 ) 


In fact, it is easy to show that all spectral forms depending 
on a single Integral scale parameter L must take the form 
of Eq. (3.17). 

When Eq. (3.17) is substituted into Eq. (3.15)5 we have 


£n[p(Sj ,825. — ^ 


N 


- Un[(a^L)^'' Pj(L) (L) . . .F (L) ] 


a^L 


Sj S2 

_P,(L) F2(L) 


+ . . . + 


Fj,(L) 


or 
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£n[p(Sj = 


|n £n(a^L) + In Fj(L) + ^,n +...+ £n Fjg(L) 


+ 


1 

a^L 



F^(L) 



F^CL) 


1 ) 


( 3 . 22 ) 


Equation (3.22) Is to be maximized with respect to o and L. 
However, from the functional form of the right-hand side 
of Eq. (3.22) we see that it Is more convenient to treat the 
right-hand side of Eq. (3.22) as a function of the two para- 
meters (a^L) and L instead of a and L. For a given observa- 
tion (l.e., sample) of our perlodogram ( S i , S 2 , . . . , Sj^j) , the 
values of (a^L) and L that maximize the right-hand side of 
Eq. (3-22) also are the values of a and L that maximize it ~ 
see, for example, p. 257 of Freeman [45]. Let us therefore 
differentiate Eq. (3.22) with respect to (a^L) and L and 
set the resulting expressions equal to zero: 


9£n[p(S^ ,82 , . . . ,Sj^; a^L,L)] 
9 (o^L) 


N 


1 


a^L (a^L)^ 


Si 


+ 


+ . . . + 


Fj(L) F2(L) 


F^(L) 


-N 


N 


(a^L)^ 


^ I 


S. 

J 


N F.(L) 

J=1 J 


1 N 

I 

(ci^L)^ j=l 


s. 

J fr 2 T 

FTiry ° ^ 

J 


( 3 . 23 ) 


and 
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ailnCp (Sj ,3^ 5 . . . a^LjD] 

3L 


^ dF^(L) ^ dP^CL) 


+ 


+ . . , + 


1 dF^(L) 


P,(L) dL PjCL) dL 


P„(L) dL 


a^L 


dP,(L) dP,(L) dP^^CL) 

+ +. . . + 


P^(L) dL 


P?(L) dL 


P^(L) dL J 


dP. (L) 

■ -J- -- 


j=l P (L) dL 

J 


S . 

sL 


dP. (L) 
J L 


N 

I - 

j=l a^L F?(L) dL 

J 


1 N 

— I 

CT^L j=l 


^ to Pj(L) 


F/L) 

0 


- a^L 


(3.24) 


where parametric dependence of p ( S ^ , S 2 , • . . , on a^L and L 

has been Indicated in the arguments in the left-hand sides 
of Eqs. (3.23) and (3-24) for clarity. Setting Eq . (3-23) 
equal to zero yields our likelihood equation for a^L: 


L = 


1 ^ 

- y 

N ^ 


S. 


J=1 F.(L) 

tJ 


(3.25) 


Setting Eq. (3*24) equal to zero yields our likelihood equa- 
tion for L after substitution of Eq. (3.25): 


N 


I 

j = l 


d 

dL 


£n P. (L) 

lJ 


F.(L) 

J 


1 N 

- y 

N .4 

1=1 



F^(L) 


(3.26) 
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Equations (Z.2S) and (3,26) are a pair of likelihood 
equations for L and o^L that involve the periodogram values 
j =• 1^2^,.. as parameters. Equation (3.26) involves 
only one unknown^ Ly and therefore may be solved first. 

Once L is obtained from Eq. (3.26) y Eq, (3.25) may be solved 
for a^. If we consider our turbulence sample w(t) to be a 
function of time of duration T, then the values of Sj are 
those obtained from Eq. (3.6) at values of f^ = j/T. 


S . 

3 


I 

T 


rT/2 

w(t) 

i 

-T/2 


^-i2i^j t/T 


dt 


2 


(3.27) 


and the values of F^(L) arey in the cases of the von Karman 
transverse and longitudinal spectray the values given by 
Eqs. (3.20) and (3.21) after substitution of 


k . = d/T . (3. 28) 

3 

In interpreting the resulting values of L as integral scales y 
care must be taken to properly account for the speed of the 
aircraft. A method for solution of Eq. (3.26) is described 
in Appendix F. 

Discussion. We can gain some insight into the general 
likelihood equations (3.25) and (3.26) by examining condi- 
tions that their solutions must satisfy. If we divide Eq . 
(3.25) by a^L and substitute Eq. (3.17) into the resulting 
expression, we obtain after minor rearrangement: 


1 

N 



E{S . } 



(3.29) 


Furthermore, if we divide Eq. (3.26) by Na^L and combine 
the resulting expression with Eq . (3.25) and define 


G. (L) = ^n F. (L) 
J dL j 


(3.30) 


we obtain 



1 

N 


S. 

E{ S , } 

L J 


1 


0 


(3.31) 


N 


■ I G (L) 


Equations (3*29) and (3.31) both are of the same general 
form — i.e., solutions L and o^L to the likelihood equations 
( 3 . 25 ) and ( 3 . 26 ) set weighted averages of [Sj/E{Sj}] - 1 
equal to zero. Since the standard deviation of the exponen- 
tial probability density of Eq . (3.10) is equal to its mean 

value E{Sj} = 1/A j , we see that the standard deviation of 
each quantity [Sj/E{Sj}] in Eqs. (3.29) and (3.31) is equal 
to the same value of unity. 


Finally^ we note that Eqs. (3*25) and (3.26) can be 
written in integral form 


a^L 




S(k) 

P(k;L) 


dk 


(3.32) 


|.k 


2 


_d_ 

dL 


Zn F(k;L) 


S(k) 

F(k;L) 




S(k’ ) 
F(k';L) 


dk' 


dk = 0 , 

(3.33) 


where continuous wavenumber k has taken the place of the 
discrete index j, and where S(k) = Sj and P(k;L) = Pj (L) . 


Likelihood Equations for von Karman Transverse 
and Longitudinal Spectra 

To specialize Eqs. (3.25) and (3.26) to the von Karman 
transverse and longitudinal spectral forms, we require the 
functions Fj (L) defined by Eq. (3.17) and the derivatives 
of their logarithms with respect to L, Eq . (3.30): 


von Karman transverse 

1+188.75 L^k! 

P,(L) s F(k,;L) = (3. S'!) 

^ ^ [l+70.78L^kn“''^ 

J 
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G (L) = ^ to F(kj;L) 


, 117.97L^k2 (1-188. 75 L^k?) 

-L J J 

^ (l+ro.78L2k2)(i+i88.75L2k?) 

J J 


von Karman tongitudinat 


F,(L) E F(k,;L) = ^ 

^ [l+ 70 . 78 L^k?]^^® 

J 


G.(L) E ^ .n P(k.;L) = - i 


117.97 L^k" 

V 

1+70.78 L^kj 


The spectra, Eq . ( 3 • 17 )j associated with the Fj(L) above are 

two-sided spectra satisfying 


- o^L F(k;L) dk , 


where k is wavenumber in cycles/unit distance, and L is the 
integral scale of the longitudinal spectrum which is twice 
the integral scale of the transverse spectrum — e.g., p. 
425 of Houbolt [^^]. The exact values of the constants in 
Eqs . (3.34) to ( 3 - 37 ) are the left-hand sides of 


25 it - = 70.78 

r(ii/ 6 ) 


200^ r(V3) . ^ 

^ r(ii/ 6 ) 


if^TT = 117.97 

^ [r(ii/ 6 ) 

— see p. 83 of Ref. I8. 
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D'isouss'ion, The von Karman transverse and longitudinal 
functions F(k;L) of Eqs. (3*34) and (3.36) are shown in Figs. 
l8 and 19, and the weighting functions Gj(L) (multiplied 
by L) of Eqs. ( 3 * 35 ) and ( 3 * 37 ) are shown in Figs. 20 and 21. 
We note first from Pigs. l8 and 19 that the functions 
F(kj;L) = P-j (L) appearing in Eq. (3.26) are Independent of 
L in the neighborhood near k = 0. Hence , spectrum samples 
Sj near k = 0 contain no useful Information for obtaining L 
from the likelihood equation (3.26) — this fact is reflected 
in the weight function Gi (L) in that Gj(L) is zero near 
k = 0. We further see tnat 


G. (L) 
0 


d 

dL 


In F . (L) 
J 


1 

F. (L) 


dF . ( L) 
J 

dL 


(3.^2) 


changes sign at kL = (l88.75)~^ for the transverse case shown 
in Fig. 20. However, from Eq. (3-^2), we see that Gn- (L) is 
zero at values of L where dF.?/dL is zero; that is, at points 
where Fj (L) is independent or L. Here again, we see that 
Gj(L) provides zero weighting to the values of k in Eq . 

( 3 . 26 ) where Fj (L) is independent of L. Finally, we note 
from Pigs. 20 and 21 that for large values of k, values of 
Gj(L) approach a constant value. Correspondingly, in Pigs, 
lo and 19 , we see for these same large values of k that 
variations in log Fj(L) with L are Independent of k. 

Finally, there is nothing in the behavior of the likelihood 
equations or the functions shown In Figs. I 8 to 21 to 
Indicate that solutions obtained for L and a^L should be 
particularly sensitive to the hypothesis that values of Sj 
are mutually Independent; hence, we would expect the 
likelihood equations (3*25) and ( 3 . 26 ) to also yield good 
results for nonstationary records with slowly varying non- 
statlonary behavior 


The statistical confidence of estimates of L and cr^L 
will be discussed in a later section of this report. 









LGi(D 



kL 

FIG. 20. FUNCTION LGj(L) GIVEN BY EQ. (3.35) FOR von KARMAN 
TRANSVERSE SPECTRUM. 
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CONSTRAINED LEAST-SQUARES ESTIMATION OF TURBULENCE 
AUTOCORRELATION FUNCTION PARAMETERS 


General Approach 

In the method for estimation of the Integral scale and 
intensity described in Sec. 3-^3 it is assumed that the tur- 
bulence time histories are drawn from a stationary, Gaussian 
random process possessing either the von Karman transverse 
or longitudinal power spectral forms. As we have discussed 
in Sec. 1 of this report, many turbulence histories recorded 
in practice have an additive low frequency (long wavenumber) 
component superimposed on what may be described as ordinary 
von Karman turbulence — e.g., see Pigs. ^ through 6. Con- 
sequently, the likelihood equations derived in Sec. 3 should 
not be used with such records. 

When a nonnegllglble fraction of the mean-square velo- 
city of a record is contained in the low frequency component 
Ws(t), we have decided to use. Jointly, the autocorrelation 
function and the power spectral density of the record to 
estimate the Integral scale and intensity of the von Karman 
component. After examination of the autocorrelation func- 
tions computed from a number of velocity histories recorded 
in the MAT Project [3^7], it became evident that over the 
range of the delay variable say where the auto- 

correlation function of the von Karman component Wf(t) is 
nonnegllglble, the autocorrelation function of the low 
frequency component Wg(t) could almost always be represented 
with reasonable accuracy by a low-order polynomial — e.g., 
see Fig. 5 of Ref. 19- A reliable model of the power 
spectral density of the low frequency component Ws(t) of 
comparable generality and simplicity is not immediately 
evident. Therefore, let us express our autocorrelation 
model as 

. m 

<K5) ^ CT?. 4>„(5;L) + I a.e^ , , ('t-D 

1=0 ^ 

where Of = E{a|-} is the mean-square value of the von Karman 
"fast" component Wf(t), is the appropriate (trans- 

verse or longitudinal) von Karman autocorrelation function 
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normalized so that 0j^(O;L) = 1, L is the Integral scale of 
the von Karman component, and 


0 


(O 


m 

= I 

1=0 


a-C' 




(4.2) 


Is the polynomial approximation to the autocorrelation func- 
tion of the slow component Ws(t) which Is valid over 
and which can contain odd as well as even powers of 
Equations ('4.1) and (4.2) are consistent with Eqs. (I.l8) 
(1.29) discussed earlier. 

We shall evaluate ^ and L by minimizing the Integral 
squared-difference E between the autocorrelation function 
model, Eq. (4.1), and the empirically obtained autocorrelation 
function R(C) that Eq. (4.1) represents: 


E 


A 





m 


1=0 


a . f ■ 
1 


d^ 


(4.3) 


where the minimization procedure will be constTained as 
follows. By definition, the "slow" turbulence component 
Ws(t) contains predominantly low frequencies In comparison 
with the von Karman "fast" component wf(t). This comment 
suggests that there usually will exist a wavenumber k£ such 
that for values of k^k£, the wavenumber spectrum of a tur- 
bulence record w(t) will be, for practical purposes, dominated 
by contributions from the von Karman component wf(t) only. 

In this wavenumber region, we therefore may use the likeli- 
hood equation (3.25) which we now write as 


L 


1 N 
— y 

N ^ 


s^. 


= 0 


j=l Fj(L) 


(4.4) 


to constrain the minimization of E. In using this constraint , 
we shall consider Eq. (4.4) to determine L as a function of 


a^. Furthermore, only the spectrum values for which 

kj>_ki will be used in Eq.(4.4); that is, the summation over 
3 in Eq. (4.4) will include only the wavenumbers k^^k^. 
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We may explicitly include the constraint, Eq. (4.4), in 
our formulation by rewriting Eq. (4.3) as 


E = 


'H 


m 


R(C) - <^^[?;L(o^)] - I a-C^i d^ , 


K 


(4.5) 


1=0 


where, for simplicity, we have used the abbreviated notation 


a 


2 



(4.6) 


and where, now, L is considered as the function of deter- 
mined by Eq. (4.4) from wavenumber spectrum values Sj for 
kj^k£. Trade-offs between choices of and m are discussed 
in Appendix G. A method that extends the model of Eq. (4.1) 
to the entire interval ’ thereby permitting computa- 

tion of the Integral scale and power spectrum of the slow 
component Ws(t), is discussed in Appendix H. 

Discussion. By this juncture, the reader has undoubtedly 
asked why the second likelihood equation (3-26) is not being 
used for the range of values kj_>kj^ to determine L. Equa- 
tions ( 3 - 25 ) and (3*26) both assume that the spectrum values 
Sj are obtained from turbulence sample functions possessing 
von Karman spectra. Thus, in choosing k^, we must be 
reasonably sure that for all kj>^kj^, the contribution from 
Ws(t) to this portion of the spectrum is negligible. By 
careful examination of spectra measured in the MAT Project 
[55], we have concluded that for a large fraction of atmo- 
spheric turbulence records, such values of k£ fall at the 
approximate location of the knee of the von Karman portion 
of the spectrum. Let us examine the behavior of the likeli- 
hood equation (3.26) for cases where the smallest value of 
j in the summation occurs slightly above this value at a 
wavenumber kj « 1/L. Examination of Figs. 20 and 21 shows 
that for values of k^l/L, Gj(L) is approximately equal to 
the constant value 

G,(L) - - ^ , k,>L"* , (4.7) 

J 0-Lj j — 

where the value of -5/(3L) can be verified for both the 
von Karman transverse and longitudinal spectra from Eqs. 

(3*35) and (3.37). If we divide Eq. (3.26) by N and sub- 
stitute the above asymptotic value of Gj(L) into the re- 
sulting expression, we obtain 
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N 


“ 3L { N ^ 


j=l 


J 

Lp.(l) 

J 


1 

N 


N S. 

I ^ 

1=1 F.(L)J 


= 0 


or 


_1 

3L 


i .1 

- J=1 


S. 


Fj(L) 


1 “ 


N S. 

1 Y 1 

N ^ 


1=1 F. (L) 
1 


= 0 


(4.8) 


which is salisfied identically for any value of L. Hence, 
when the smallest usable value of kj Is L~^ or larger, 
the likelihood equation (3.26) Is satisfied identically for 
all values of L and therefore is useless. For somewhat 
smaller values of k£, Eq . (3-26) will still yield unreliable 

results. Hence, in essence, we have replaced Eq. (3.26) by 
minimization of the parameter E described by Eq . (4.5). 

On the other hand, the first likelihood equation (3-25) 
is perfectly well behaved when the only usable values of kj 
are the values kj>L In fact, for the records that we 

shall discuss in the applications portion of this work, the 
likelihood equation ( 3 . 25 ) yields results with excellent 
statistical reliability. It therefore has been retained as 
the constraint L = L(a^) in the minimization of Eq . (4.5). 

There exist two additional reasons for using Eq. (4.4) 
as a constraint in minimization of the Integral squared- 
difference E given by Eq . (4.5). First, we note that when 

the degree m of the polynomial of Eq. (4.2) is taken too 

large, the von Karman autocorrelation function 0j^(^;L) 
can be represented quite well by the polynomial of Eq . (4.2). 
When this is the case, minimization of F — given by Eq. 

(4.3) — leads to set of ill-conditioned equations for L, 

crp, and ao jai , . . . ,aj^. On the other hand, when Eq. (4.4) is 
Included as a constraint in the minimization of E — as 
indicated by Eq. (4.5) ~ the von Karman component 

^ ^ remains dependent on only one parameter 

In the minimization; consequently, the range of shapes that 
^ L ( ) ] can take on in the minimization is greatly 
reduced, and for a given value of m, the minimization be- 
comes much better conditioned. This is a very important 
consideration when working with empirical data such as the 
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empirical autocorrelation function R(^) In Eq. (^.5). We 
shall see In the applications portion of this work how this 
Improvement in conditioning yields improved results In 
certain specific cases. 

The second reason for using the constraint L = L(a^) 
of Eq. (4.4) is that we shall want to estimate the wavenumber 
spectrum of the ’’slow" component WgCt) by subtracting our 
best estimate of the spectrum of the von Karman component 
from the overall spectrum of w(t). Since the equation of 
constraint (4.4) was obtained directly from the wavenumber 
(or frequency) domain and represents an optimum estimate 
obtained from the von Karman portion of the spectrum, we 
would expect to obtain better results for the spectrum of 
Wg(t) when Eq. (4.4) Is Included as a constraint In the 
minimization of E. 

Finally, we should comment on why it Is possible for 
the minimization of F — as described by Eq. (4.5) and the 
constraint equation (4.4) — to yield potentially better 
estimates of and L than use of the pair of likelihood 
equations ( 3 . 25 ) and (3.26) over the range of values of 
kj^k£ where the wavenumber spectrum Is dominated by the 
spectrum of the von Karman component wf(t). In setting up 
the quantity E represented by Eq. (4.5) to be minimized, 
we have included the representation, Eq . (4.2), of the auto- 
correlation function of the low-frequency turbulence com- 
ponent, Wg(t). Equivalently ^ we have included a vepvesenta- 
tion of some of the information about the ''slow” turbulence 
component Wgit) that would be found in the wavenumber 
domain in the region k<ki. We also have Included Information 
about the von Karman component that would be found In the 
same low-wavenumber region k<k^. The reason we have chosen 
to work with the autocorrelation function representation 
of Eq. (4.1) rather than with a comparable representation 
of the spectrum of w(t) = Ws(t) + Wf(t) Is that It Is a 
simple matter to make intelligent a priori choices of 
and m for our representation of 4>w (?) given by Eq. (4.2) 

^ 3 

which should lead to good results when the Integral squared- 
difference given by Eq. (4.5) Is minimized. A representation 
of the wavenumber spectrum of Wg(t) of comparable generality, 
simplicity and amenableness to analysis Is not obvious. 
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Derivation of Algebraic Equations for 
Autocorrelation Function Parameters 

Let us now derive the equations whose solution yields 
a minimum value of E. That is, from Eq. (4.5), we wish 
to derive a set of algebraic equations whose solution for 
and ao ,ai 5 . . . ,aj^ minimizes the value of E given by Eq. (4.5), 

where Eq. (4.4) with e serves to define L as a function 
of for use in the resulting set of equations. These 
equations are obtained by setting the derivative of Eq. (4.5) 
with respect to ao ,ai , . . . ,ain and equal to zero, where L 
is considered as a function of : 


"bE 

8a 


= 0 , j=0,l,...,m 


bE 


= 0 . 


(4. 9a,b) 


9a" 

Differentiating Eq. (4.5) with respect to aj , we obtain 


bE 

8a. 

J 


= -2 


’H 


m 


R(5) - $ [C;L(a")] - I a } 5^ <35 , 

1=0 ‘ 


j = 0 , 1 , . . . , m 


(4.10) 


Furthermore, differentiating Eq. (4.5) with respect to 
while treating L as a function of we obtain 


bE 

ba^ 



;L(a^ ) ] 



X 



8L 


dL 

da^ 


+ (l)j^[^;L(a")] 


dC . 


(4.11) 


Setting Eqs. (4.10) and (4.11) equal to zero yields the set 
of m + 2 nonlinear algebraic equations for a o , a i , . . . ,a^^ 
and which can be written as 
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£j(tj^C5;L(0^)] ds + J d5 = 5^R(5) dj. 


= 0 ^ 1 


J=0,l,...,m (4.12) 


and 


) 2 dL , , r. -r/ 

+ (t)^[^;L(a^ ) ] 


3L 


da 


X (). j ^[ 5 ; L ( a 2 )] ds 


m 

+ I a, 

1=0 


) ^2 dL 


3L 


da' 


+ <t>jj[5;L(o^)]}5^ dC 


■H I 2 ^ ^ 


a' 


8L 


dL 

da^ 


+ (|)j^[C]L(a2)] } R(^) d^ , (4.13) 


where, from Eqs. (4.4) and (4.6), L(a^) is defined by the 
equation 


-, N S. 

cj 2 = 1 y — i 

^ M Z. 


” j=l LF (L) 

t) 


(4.14) 


where, as Indicated earlier, only values of Sj dominated by 
the von Karman component are Included within uhe summation 
In Eq. (4.14). 

For the von Karman transverse and longitudinal spectra 
of Eqs. (3-18) and (3.19), let us define the normalized 

spectrum 0(k), where 

k = Lk , (4.15) 
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by 


= a^L 


(4.16) 


where the overbars denote normalized spectra and normalized 
wavenumber. Then, for either transverse or longitudinal 
spectra, the corresponding transverse or longitudinal auto- 
correlation function is 




$^(k) dk 

J ^ 


a^L I $^(Lk) dk 

_oo 


= 4-g(k) e^2rt5/L 


(^.17) 


If we define a normalized length measure by 

5 - ?/L ('t. 18 ) 

and a normalized autocorrelation function by 




(4.19) 


then, according to Eqs. (4.17) to (4.19)s we may express 
4 >k(^) as 


(j>j^[^;L(a2)] E (J)j^(^) = (f)j^(?/L) = cj)j^(0 


(4.20) 
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Consider, now, the derivatives of <j)j^ with 
L which appear In Eq. (4.13). According to Eq. 
have 


3L L" 


hence, from Eq. (4.21), It follows that 

04)j^[C;L(a2)] af 

dL dC 9L 

= _ _L 

df 


Turning to the integrals that appear In Eq. (4, 
define 


X^(L) ^ 



J 

0 


— (})^[C;L(aM] dC 

9L ^ 


('CtJ'/I-J 

j 4 (l)j^(0 (J)j^(C) d? 


where Eqs, (4.l8), (4.20), and (4.22) were used 
the second line, and where 




d<j>j^(?) 

d¥ 


The second Integral that we require Is 


respect to 
(4,l8), we 

(4.21) 

(4.22) 
13)j let us 

(4.23) 

In going to 

(4.24) 


100 



0 




& 1 

L 


•H 




J “ ¥|(f) df , 


(i).25) 


where Eqs. (4.18) and (4.20) have heen used in going to the 
second line. Next, we require 


JaCj.L) ^ 


i 3L 


dC 


^ 4)^(0 dC , j = 0,l 


,m 


(H.26) 


where Eqs. (^.l8) and (^.22) have been used. Prom Eq. (4.13), 
we also require 


X,(J,L) 


4 


1 



<l)j^[5;L(a^)] d5 


J <(>k(C) dC , J = 0,l,...,m , (4.27) 


where Eqs. (4.18) and (4.20) have been used. We next 
require 
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ZJL) ^ I “ — ^ R(5) dC 


R(0) 


3L 


R(0)L^ i 


S <tK(?/L) R(C) d5 


1 


R(o) i 


C„/L 

? .(.g(C) R(L5) dc , 


Where Eqs. (4.18) and (4.22) have been used. We 
require 


T,(L) = 


A 1 


R(0)L 


‘H 


c()j^[^;L(g")] R(0 dC 


'H 


R(0)L 0 


(|)j^(C/L) R(C) dC 


1 


R(0) ! 


^rr/L _ _ 

(p^(^) R(L^) dC , 


where Eqs. (4.18) and (4.20) have been used. Two 
integrals are required for use in Eq. (4.12): 

d5 


1+J 


J = 0,1, ... ,2m 


(4.28a) 


(4.28b) 

further 


(4.29a) 


(4.29b) 

additional 


(4.30) 
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and 


-r.(J,L) i 


R(0)L^'^J I 


5^ R(5) d? 


(it-Sla) 


R(0) 


5h/^ -1 - - 

R(L5) d5 , J=0,l,...,m (4.31b) 


where Eq. (4.l8) has been used. The above quantities 
through Jg are dimensionless. 

Using the definitions of Eqs. (4.23) to (4.31)a we can 
rewrite Eqs. (4.13) and (4.12) as 


a 


2 



dL 

da^ 


J^(L) + L J^(L) 


+ 


m 


I 

1=0 



da^ 


■ 


= -R(0) J,(L) + R(0)L I AL) (4.32) 

da" ' ' 

and 

a"L IJJ,L) + I a. Jj(i+J,L) = R(0)L J 3 (J,L) , 


j=0,l,. . . ,m . (4.33) 


Equations (4.32) and (4.33) for j=0,l,...,m constitute a 
set of m + 2 nonlinear simultaneous algebraic equations for 
a" and ao , ai , . . . , a^^. For a given value of a", Eq. (4.l4) 
(inverted) determines the quantity L(a"). Equation (4.l4) 
also determines dL/da" : 
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L 



5 


(4.34) 


dL 

da^ 



1 N 

- y 

N .4_ 
J=1 


[F. (L)+L 
J 


FI (L)]S. 
.1 >1 


- 1 


L^F? (L) 
J 


which yields dL/dc^ for a given value of L(a^). 

Matrix Form of Algebraic Equations for 
Autocorrelation Function Parameters 

For numerical solution, the set of m + 2 equations 
(4.32) and (4.33) can be written In the matrix form 

AiiCaOyi + +■••+ 

A2i(a^)y;^ + AgjCaOyg +•••+ 

+ * m + 2 , 2 (<^') y 2 +•••+ \+2 , m +2 ' ) ym +2 = %+2 

(4. 

or, more concisely, as 
m+2 

^ k=l,2, . . . ,m+2 (4. 

where 


^S. ^ ^1-2 1=2,3, ... ,m+2 


5 R(0) 


-O' 


dL 

da^ 


J,(L) + L J ,(L) 


x^ = R(0)L Jg(k-2,L) for k=2 , 3 , . . • ,m+2 


(4. 


(4. 


(a^) 

35a) 

35b) 

36) 

37) 
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A = -o^ J,(L) + L J (L) 

A,. H -o^ J,(Jl-2,L) + lit-2,1,) , 

da^ 

a=2 5 3 s . • • ,ni+2 

E L I-Jk-2,L) , k=2,3, . . . ,m+2 

= L^”l r^(k+Z-H,L) , k=2,3, . . . ,m+2; «=2 , 3 , ■ • • ,m+2 

(It. 38) 


Method of Solution 


Scaling of matrix equations. The m + 2 equations (4.35 
were scaled before their solution was obtained — using the 
method suggested on pp . Il8— 119 of Hamming \^44~\. This 
scaling Is carried out as follows: 


1. Consider a new array of "coefficients" which is the 
array In Eq. (4.35) with the right-hand side added as 

an additional column 


11 

A 

^12 

^l,m+2 


'21 

• 

A 

22 

• 

• 

^2,m+2 

• • 

• 

^2 

• 

• 

m+2 5 1 

• 

^m+2,2 

• 

■^m+2 ,m+2 

• 

^m+2 


2. Call the above new array where 1 ranges from 

1 to m + 2 and j ranges from 1 to m + 3* 

3. Form a new array Aj_^ from the array A - where for 
every 1=1 , 2 , . . . ,m+2 and j =1 , 2 , . . . ,m+3 



r . +c . +M 
2 ^ J 


A I . 
IJ 


('j.sg) 
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Formulas for r^, Cj , and M will be given In Item 6 below. 

Let yi>y23 • • • aym+2 solution to a new set of 

matrix equations which may be written in terms of our ori- 
ginal notation as 


2 


r . +c . +M 
1 J 



r.+c ,o+M 
, 1 m+3 



(4.40) 


M 

5. Multiplication of every element by 2 does nothing 

r* . 

to the solution. Multiplication of every row by 2 ^ does 
nothing to the solution. The effects of multiplication of 

the rle:ht-hand side by 2 can be eliminated by dividing 

Cm+ 3 

every element by 2 . Therefore, the relationship between 

yj and y^ is 


c 

2 


3 


2^m+3 


y. 


y 


(4.4la) 


In other words, for every J the solution y.j of 
equations (4.35) is obtained from the solution 
(4.40) by 


the original 
yj of Eqs. 


c — c 

yj = 2 ^ y^ , j = l,2, . . . ,m+2 . (4.4lb) 

6. To determine values for rj_, cj , and M we first form 
for every pair of values l,j where i=l , 2 , . . . ,m+2 and 
j=l,2, . . . ,m+3 the quantity 


logi„2 


(4.42) 


where |Aj^j| denotes the absolute value of Aj_ . , and where 

3 

Ai item 1. Values of M, r^, and c^ are then 

computed by 
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M 


-1 


(4.43) 


(m+2)(m+3) 


m+3 

I 


j=i 


m+2 


I 

1=1 


d 


IJ 



m+3 



^1 ~ m+3 

I 

J=1 

(d. .+M) 
ij 



m+2 



1 

m+2 

1 

1=1 

(d. .+M) 

(4.45) 


Initialization of computation. Steps in preparing for 
solution of the matrix equations (^.40) are: 

1. Compute the unsmoothed power spectral density of the 
record of Interest. Call the (unsmoothed) computed spectrum 
values Sj_. Pick lower and upper wavenumbers k£ and which 
define a wavenumber region k^<k<k^ where we are confident 
that the spectrum is dominated by the von Karman component 
of the turbulence. For values of Sj[ that fall within this 
wavenumber region, determine for a set of equally spaced 
values of Lj the quantity 


2 - 


= a 


(Lj) 


1 ^ 

= ~ y 

^ 1^1 


^1 

F. (L. ) 
1 J 


(4.46) 


as defined by Eq. (4.l4), where Fi(L) is given by Eq. (3-20) 
for transverse (vertical or lateral) records and by Eq. 
(3-21) for longitudinal records. When a^(Lj) is computed 
for a range of equally spaced values of Lj , we have as a 
numerical function of the integral scale Lj . This tabula- 
tion is to be considered as determining L as a function of 
a^. Values of L falling between tabulated values of Lj 
are obtained by Interpolation. 


2. Values of dL/da^ also are required as a function of 
. In our computations, we used the approximation 


dL 

da^ 



L. 

_1 

a? 

3 


(4.47) 
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where the values of L and In the right-hand side of Eq. 
(4.47) are those computed by Eq. (4.46). A more accurate 
procedure would be to compute dL/da^ by Eq . (4.34) for the 
same values of Lj used in Eq. (4.46). 

3. From the autocorrelation function R(^) of the 
record under consideration, and for the same values of Lj 
used in Eqs. (4.46) and (4.47), the Integrals of through 
Jq given by Eqs. (4.23) through (4.31) are computed numer- 
ically. Each of these Integrals is a function of the 
integral scale L. Thus, Ii through le is to be computed for 
each of the above mentioned equl-spaced value of Lj . The 
values of the Integrals are then tabulated. Computation of 

these integrals requires the functions cj)|^( 5 ) and 4(5) which, 
for von Karman transverse and longitudinal spectra, are: 

for von Karman transverse spectra (vertical and lateral 
components ) : 


<1>k(5) 

4(5) 


2/3 


r(i/3) 


( 65 ) 


K^3(65) - I 5 (65) 


3 


(36)"^' K^^3(30 - I (3^)''^ K_,/3(30 


8 /" o ir ^ 1/3 


2^^^r(i/3) 


(4.48) 


(4.49) 


for von 

Karman 

longitudinal sp 

4 ( 5 ) = 

22/3 

(65)^^' 

^ 3 ( 65 ) 

r(i/3) 

4(5> = 

_ 22 /Sg 

(65)*^' 

K_2/3(55) 

r(i/3) 

where 




6 

A 2/? r(ll/ 6 ) 
5 r(4/3) 



ctra : 

(4.50) 

, (4.51) 

(4.52) 
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and ^_ 2 /^ modified Bessel functions of the second kind 

of order 1/3 and -2/3, and f(*) is the gamma function. 

Values of the Bessel functions were read off of the tabula- 
tion on p. 228 of Ref. 45, where we note that K_^^(x) = 

K 2 / 3 (x). Values of ^ were tabulated rather than values 

of because has a singularity af ^ = 0 whereas 

? tends to zero as 1" 0. Values of Ji through Jq 

required for L / Lj were obtained by Interpolation between 
the tabluated values of the T’s. 

Solution for and ao ^ ai ^ ^ a^. Once L(a^), dL/da^, 

and 1 1 through Ib are tabulated for a predetermined set of 
equl-spaced values of L = Lj , the solution to the set of 
equations (4.40) can be obtained. This solution yields the 
solution to the set (4.35) by Eq. (4.4lb). Although the 
scaled equations (4.4o) were actually the set solved, we 
shall describe the solution procedure in terms of Eqs. (4.35). 

The coefficients and the right-hand sides Xj^ of 

Eqs. (4.35) are functions of the unknown variance = yi 
of the von Karman component. Thus, Eqs. (4.35) were solved 
by trial and error. Prom a plot of the autocorrelation 
function of the record R(^), a rough estimate of 
is easily obtained. This value of determines L = L(a^ ) 
and dL/da^ which are obtained by Interpolation from pre- 
determined tabulated values of these quantities. Ii through 
Tb also are determined by Interpolation from predetermined 
values tabulated as functions of L = Lj . Consequently, 
once the initial value of has been chosen, the coef- 
ficients AjL j, i , J = 1 , 2 , . . . ,m+2 and right-hand side x^_, 

1=1 , 2 , . . . ,m+2 in Eqs. (4.35) are determined. Equations 
(4.35) are then solved, the solution being a new value of 
and ao ,ai , . . . ,ajji [see Eq . (4.36)]. Using the new value of 
, the coefficients Aj_j , i , j = 1 , 2 , . . . ,m+2 and right-hand 
side xi, 1=1 , 2 , . . . ,m+2 are again evaluated, and the set of 
Eqs. (4.35) is solved again. The new solution yields a 
new value of and values of ao ,ai , . . . ,am. By comparing 
the solution values of with the input values of for 
these two solutions, a new trial value of is chosen and 
the set of Eqs. (4.35) is solved again after evaluation of 
the coefficients and the right-hand side using the new 
trial value of c^. In carrying out this procedure, we 
terminated the process when the input and solution values of 
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were in agreement to three significant figures. The final 

values of e q|. and L(a^) determine the spectrum of the 
von Karman component of the turbulence. 
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VARIANCE OF MAXIMUM LIKELIHOOD ESTIMATES OF von KARMAN 

TURBULENCE PARAMETERS 


In Sec. 3, a pair of equations (3*25) and (3.26) was derived 
whose solution yields maximum likelihood estimates of o^L and L 
from spectrum samples Sj , j = 1^2, These spectrum samples 
are to be computed from turbulence velocity records whose 
(expected) power spectral densities are assumed to be of the 
von Karman transverse or longitudinal forms . Since Individual 
turbulence records are stochastic functions of time^ estimates 
of a^L and L computed from such records also are stochastic. 

Thus, in any particular application the solution cr^L, L to the 
pair of equations ( 3 * 25 ) and (3.26) Is a two-dimensional random 
variable whose joint probability density function depends on the 
duration of the turbulence sample function from which the estimate 
of a^L and L Is obtained. In this section, the probability 
density functions of maximum likelihood estimates of o^L and L 
are discussed. In addition, explicit formulas are given for the 
squares of the ooeff-io'ients of vaviat'ion — i.e., relative vari- 
ances — of a^L, L, and a^. 


Asymptotic Forms of Probability Density Functions 
of Turbulence Parameter Estimates 

In the "standard" class of problems that employs maximum 
likelihood estimation: each of the Individual probability density 

functions A^-e , j = 1,2, •••,N In Eq . (3.11) Is identical. 

However, this Is not true in the present application because the 
values of the parameters Aj depend on J — as may be seen from 
Eqs. (3.12), ( 3 . 17)5 ( 3 . 20 ;, and (3-21). Hence, care must be 

exercised In applying the "standard" results to the present 
problem . 

The "standard" result that Is our main Interest Is the 
asymptotic (large sample) form of the joint probability density 
function of our estimates of a^L and L obtained by solving Eqs. 
( 3 . 25 ) and ( 3 . 26 ). In the standard class of problems mentioned 
above, such maximum likelihood estimates are jointly normally 
distributed In the large sample limit — see, e.g., p. 55 of 
of Ref. 38 or p. 155 of Ref. 46. In the present case, we shall 
appeal to the form of the multidimensional large sample result 
cited on p. 155 of Ref. 46, which avoids the troublesome problem 
of bias, and which is the multidimensional extension of the 
approach used by Cramer on pp . 550 to 504 of Ref. 8. 
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To Justify use of the large sample asymptotic normal ■ form 
in the present "nonstandard” applicatlona we must show as the 
total duration T of our turbulence sample function approaches 
infinity that (1) 9£np/9(a^L) and 9£np/9L are asymptotically 
Jointly normally distributed with zero expected values at the 
true values of a^L and L and (2) sample values of 9^£np/9^ (a^L) , 
9^£np/9 (a^L) 9Lj and 9 ^£np/9L^converge to their respective 
expected values in the sense that fractional deviations from 
their expected values vanish with probability one as T ^ — see 

pp. 154 and 155 of Ref. H6 and pp . ^3 and 55 of Ref. 38. 


Consider requirement (1) first. According to Eqs. (3*23) 
and (3. 24) 5 9£np/9(a^L) and 9Anp/9L both are linear combinations 
of the Independent random variables (Sj/P-j), J = 1,2, From 

Eqs. (3.10), (3.12), and (3.17), it follows that each random 
variable (Sj/Pj) is governed by an exponential probability 
density function with the same (unknown) expected value of a^L. 
From this fact, it follows that at the true values of a^L and L 
we have E{ 9£np/9 (a^L) } = 0 and E{9£np/9L} =0. To show that the 
two-dimensional random variable { 9^/np/9 (a^L) , 9£np/9L} is 
asymptotically normally distributed as T ^ employ the two- 

dimensional central limit theorem — e.g., pp . 285-287 of Ref. 8. 
Since, in the case of Eq . (3.24), we are dealing with a weighted 
sum of Independent random variables, where the weighting function 


is 


d 

dL 


£nP . (L ) , 
J 


we first consider the behavior of the sums in 


Eqs. (3.23) and (3*24) over a typical limited fixed frequency 
range, say Af, as T ^ According to Eq. (3. 13)5 if the fre- 

quency range Af is fixed, as T ^ the number of samples within 
Af Increases indefinitely since the frequency difference between 
adjacent samples of 1/T. This behavior is true for every fre- 
quency interval Af over which the histogram S(f) — Eq . (3.6) — 


is computed. Since the slope of the weight function ^ £nPj(L) 
is everywhere finite - see Eq . (3.30) and Pigs. 20 and 21 - we 

may choose Af small enough so that JlnP.- (L) is essentially 


constant within Af. We may now apply the central limit theorem 
to the contributions of the sums in Eqs. (3-23) and (3.24) within 
each such frequency Interval Af to show that the contributions 
from each such Af are asymptotically Jointly normally distributed 
as T The sums in Eqs. (3.23) and (3.24) over all such dis- 

joint intervals then represent the sums of ( nonldent Ically 
distributed) Independent two-dimensional normal variables which 
are known to be normally distributed — e.g., pp . 212 and 316 of 
Ref. 8. Hence, as T the two-dimensional random variable 

{ 9inp/9 (a^L) , 9^np/9L} is asymptotically normally distributed 
with zero mean value at the true values of a^L and L. 
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Let us now consider requirement (2) above. Forming the 
partial derivative of Eq. (3*23) with respect to a^L gives 

a2£nrp(Sj • • -,S^;a^L,L)] 


a(a^L)^ 


N 


2 N 

^ — I 

(a^L)2 (o^L)^ J=1 


P. (L) 

J 


- a^L 


(5.1) 


forming the partial derivative of Eq. (3.23) with respect to L 
gives 

3Hn[p(Sj,S2,---,Sj^;a^L,L)] 

3(o^L)3L 


(a^L)^ j=l 


N dP.(L) 

I [P,(L)]-^^S, 


1 “ 

^ I 

(a^L)" J=1 




s . 


sTtL) 

J 


(5.2) 


and forming the partial derivative of Eq . (3.24) with respect to 
L gives 

9Hn[p(S^,S^,---,S ;a"L,L)] 


9L' 


a^L j=l 


ilnP. (L) 
dL j " " 


S. dF.(L) 
■J J 


p!(L) dL 
J 


d^ 

— £nP. (L) 
LdL^ J 


r s 




- a^L 


or 
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Let us define the fractional deviation of 8 ^ £np/ 9 ( a ^ [given by 
Eq. (5.1)] from Its expected value [given by Eq. (5.5)] by — 
l.e. , 
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A [ 8^ Jlnp/3 (g^L ) ^ ]-E{ 3^ilnp/8 (g^L) ^ } 

E{8^ilnp/8(a^L) 2} 

8^£np/8(a^L) ^ ]-E{ 8^£np/8 (a^L) ^ 

I E{9Hnp/9(a"L)"} , (5.8) 


and denote analogous definitions of the fractional deviations of 
9^ ilnp/3 (a^L) 3L and 9^S,np/9L^ by e and £3 respectively. From 
Eqs . ( 5 . 1 ) to ( 5 . 3 ) and (5*^) to (5.6), it follows directly from 
the expression for e^ given by Eq. (5-8) and the corresponding two 
expressions for and that 


= 


p 1 ^ 

— - y 

2t N 

a^L j=l 


S . 

f7l) 

ij 


- a^L 


(5.9) 


a » = 


1 N 

- y 

N . 

j=i 


§L 


s . 


PTtL)- 

J 


- a^L 


1 N . 
J=1 


(5.10) 


and 


N 


- y 

N ^ 


e , = 


J=1 




dL 


~ JlnF . (L) 
2 J 


S . 


N 


F. (L) 
: J 


- a 


I 


J=1 




(5.11) 


We must show that e^, ^3 approach zero with prob- 
ability one as our turbulence sample duration T Since the 

random variables Sj/Fg(L) are each governed by an exponential 
probability density with mean value a^L, it follows from the 
(weak) law of large numbers — e.g., p. 228 of Ref. 47 — that 
e j 0 as T -> <» with probability one. 
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To show that £2 ^ note first that, in the 

limit T the denominator of the right-hand side of Eq. (5.10) 

approaches the limiting constant value 


V 


2 


A a^L 
(N2~N^)Ak 


rN.AK 

1 

N^Ak 


d 

dL 


£nF(k;L)dk . 


( 5 . 12 ) 


This fact becomes clear when we observe (1) that the denominator 
in Eq. (5. 10) is a^L times the average value of the derivative of 
£nFj(L) with respect to L, (2) that, in the large sample limit, 
dependence of Fj (L) on the discrete variable j is replaced by 
dependence on the continuous variable k as shown in Eqs . (3.3^) 
and (3.36), and (3) that in taking the limit T ^ in applications, 
we fix the lower and upper wavenumber limits 

k = N Ak, k = N Ak (5.13) 

1 1 5 2 2 \ -J / 


where N = N2~N^ and Increase indefinitely the density of sampling 
points within the Interval so that the sum approaches the 

definite integral in Eq. (5.12) of the continuous function 


G.(L) = ^ £nF(k.;L) E ~ £nF.(L) 
J dL J ’ dL j ' 


(5.14) 


shown in Figs. 20 and 21 for von Karman transverse and longi- 
tudinal spectra. 

Consider, now, the numerator of Eq. (5.10). Here, we may 
appeal directly to the form of the (weak) law of large numbers 
given on p. 238 of Ref. 47, which applies to sequences of non- 
Identlcally distributed Independent random variables. The 
numerator in Eq. (5. 10) can be expressed as 


N 



1 ? 

^ J=i J 


(5.15) 


where the 


X . 


J 


dL 


InF. (L) 
J 


S . 

J 

fTTTT 

J 



(5.16) 
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are the random variables of Interest. Since Sj/Fj (L) Is governed 

by an exponential probability density with mean a^L (Independent 
of 0)3 It follows that the random variable { [Sj/Pj (L) ]-a^L} 
has zero mean and variance (a^L)^. Hence, the random variables 
defined by Eq. (5.16) have zero mean and variances 


a 


2 


<:3 


(a^L)^ 


dL 


^.nF . (L) 


2 


(5.17) 


A sufficient condition for to vanish with probability one as 
T 00 is — according to Eq. (5.6) on p. 238 of Ref. ^7 ~ that 


^ ^ I 

N2 ^2 J 


(g^L)^ 


N 

I 

j=l 


^ ^nF. (L) 
dL j 


( 5 . 18 ) 


vanish as T However, In this limit Eq . (5.18) approaches 


N(Nj-N,)Ak 


rN Ak 

/ 

NjAk 


d 

dL 


^nP(k;L) 


2 

dk 


(g^L) ^ 
N(k^-k^) 



JlnP(k ;L ) 



( 5 . 19 ) 


where we have substituted N = N^-N^ and have written the^ limiting 


form of the sum on the continuous function 


^ toFCkj.L) 


as 


an integral as before. The right-hand side of Eq . (5*19) is pro- 
portional to 1/N; hence, (S]^/N) is proportional in the limit to 


1//N and therefore to 1//T. therefore approaches zero with 

probability one as T -> «». The argument that must approach 
zero with probability one as T -j- °° Is carried out in exactly the 
same way as that for £ 2 . 
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The maximum likelihood est-imates of a^L and L therefore are 
asymptotically jointly normally distributed with mean values equal 
to the true parameter values (denoting estimates of a^L and L 
subscripts 1 and 2 respectively) 

m^ = m^ = L, ( 5-20 a,b) 

and with covariance matrix A whose inverse Is given by 


A 


1 1 


A 


1 2 


A 


- 1 _ 


A 


2 1 


A 


2 2 


where 


= -E 


5 ^ £np 


9(a^L)^ 


A,2 = -E 


9 ^ £n p 
9L^ 


and 


A = A 
^12 ^21 


_E I 

9(a^L) 9L) 


( 5 . 21 ) 


( 5.22 a,b ,c ) 


(p. 55 of Ref. 38 and p. 155 of Ref. 46). Expressions for the 
(negatives) of A^^, ^ 12 ’ ^21 given by Eqs . (5.5) to (5-7). 

It is easy to verify that the Inverse of A“ ^ is given by 


A 






(5.23) 


where |A is the determinant of A ^ , 


_ 1 


A 


- 1 


= A , A -A A , 
11 22 12 21^ 


(5.24) 
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and where the form of Eq. (5.23) can be verified by multiplying 
Eq. (5.21) by (5.23) and using Eq. (5.24): 


1 




iO 


0 


1 


(5.25) 


If we write the covariance matrix A in the conventional form — 
e.g., p. 295 of Ref. 8 — 


A 



pa 


2 





(5.26) 


it follows by comparison of Eqs . (5.22)^ (5.23), (5.26), and Eqs . 
(5.5) through (5.7) that we have 


I A"M 


9^£np ( 
9(g^L) " ) 


|A-'| 


and 


N 


2 

1 

v)=l 

J 



|A-' 


N 

(g^L)^ 

|A-M 


p = 


) 3^£np I 

(9(g^L) 9 l) _ 

• |A"M 

12 ' ' 


N 


I ^^nF(L) 


J=1 




N 

N I 
J=1 


dL 


(5.27) 


( 5 . 28 ) 


(5.29) 
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where from Eqs . ( 5 . 22 ), (5.24), and (5.5) through (5.7), we have 


1 


= E' 


9 ^ £np 
9 (a^L) 


E 


(a^L) 


N 

N I 

J=1 •- 


9 ^ 
d 


ZnF.(L) 
dL j 


9^£np I 
9(a2L)9L| 


N 

I 

Lj=i 


1 2 




(5.30) 


In summary , we have shown that solutions to the patr of 
ttkeZthood equat'tons (2.25) and (2.26) proV'ide estimates to the 
true values of o^L and L that as T are asymptotically 

governed by a joint normal probability density with mean values 
equal to the true parameter values a^L and L as indicated by 
Eqs. (5.20 a,b) and with covariance matrix elements given by 
Eqs. (5.27) to (5.20)^ where and of denote, respectively , the 
variances of our estimates of o^L and L, and p denotes the 
correlation coefficient of these estimates. The limiting forms 
(as T -^ «) of Oj, and p for the von Karman transverse and 

longitudinal spectral forms will be evaluated in the next section. 


Expressions for Covariance Matrix Elements for 
von Karman Transverse and Longitudinal Spectra 

Let us turn now to evaluation of the large T limiting forms 
of the sums in Eqs. (5.27) to (5.30) which apply to our asymptotic 
Joint normal distribution of the maximum likelihood estimates of 
a^L and L. We shall continue to use the notation established 
in Sec. 3 — Eq. (3.17) in particular — and at the end of this 
section we shall specialize the results to the von Karman spectral 
forms of Eqs. (3.20) and (3.21). 


First consider the expression for p given by the right-hand 
side of Eq. (5.29). If we divide both numerator and denominator 
of Eq. (5.29) by N, we shall require expressions for the two 
quantities 


t(1) a 1 
N 



d 

dL 


InF. (L) 
J 


(5.31) 


and 
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(5.32) 


(2) A 1 
N 


N 

I 

J=1 


d 

dL 


2 

inF. (L) 


Since the right-hand sides of Eqs . (5.27) to (5.30) represent the 
limiting forms of the variances and correlation coefficient of 
our estimates of a^L and L as T oo, we can, without loss in 
generality, consider only the limiting forms of the right-hand 
sides of Eqs. (5-31) and (5-32) as T Denote these limiting 

forms by . 


(1) 4 t(1) 

T->oo 


(2) A 11m j(2) 


(5-33 a,b) 


and consider I first. Here, we shall consider the record to 
be of length L where L = VT, T being the duration of the record 
and V being the speed of the aircraft used in measuring the 
record. Hence, our N wavenumber spectrum samples Sj discussed 
in Sec. 3 are spaced at wavenumber intervals of 


Ak = 1/L 


(5.34) 


[which correspond to frequency intervals of Af = 1/T — see Eq . 
(3.13)]. 

Recognizing that the right-hand side of Eq. (5.31) repre- 
sents an averag'e of ( d/dL ) ^nFj (L ) , and that as T and L ^ our 
wavenumber spacing Ak ^ 0 , we have for the limiting form of 

given by Eq. (5-33a): 


Cl) 1 a 

^ = filk J k ^nF(k;L)dk, 

NjAk 

where 

N = , 


(5.35) 


(5.36) 


and where we have used the notation of the left-hand identities 
in Eqs. (3.34) and (3-36). From Eqs. (3.35) and (3*37), we see 
that L(d/dL) 5,nP(k;L) is a function of only the product Lk; hence, 
let us define 
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H(Lk) - L ^ £nF(k;L) . 


(5.37) 


Using Eq. (5.37), can be expressed from Eq . (5-35) as 


( 1 ) _ 1 


NL^Ak 


i-N Ak 
' 2 


NjAk 


H(Lk)Ldk 


(5.38) 


NL^Ak 


l-N^LAk 


N^LAk 


H(^)d^, 


(5.39) 


where ^ = Lk . Let us define 


H ^ U” H(?). 

oo f-^CO ' ^ 


(5.^0) 


Recognizing that we are especially Interested in the case of 
"large" upper limits N^Ak and N^LAk in Eqs . (5.35) and (5.39) 
respectively, let us now rewrite Eq. (5*39) as 


( 1 ) _ 


1 

NL^Ak 


N LAk 


[H(?)-H^]d? + 


H (N,-N, )LAk 

oo 2 1 


N^LAk 


NL^A k 


1 

NL^Ak 


rN LAk 
‘ 2 


N^LAk 


H 

[H(5)-H„]d5 + , 


(5.41) 


where we have used Eq. (5.36). 

The likelihood equations (3*25) and (3-26) apply for any 
limits Nj and N 2 in Eq. (5.^1). However, generally we are 
interested in a (dimensionless) low wavenumber limit of 
kjL = NjLAk = 0, whereas the (dimensionless) high wavenumber 
limit k 2 L = N 2 LAk is usually taken to be the largest wavenumber 
for which the spectrum is still unaffected by instrumentation 
errors etc. This highest wavenumber is usually of the order of 
k^L 10, where L is the integral scale of the turbulence. From 
Eqs. ( 3 . 35 ) and (3.37), we see that 
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( 5 . 42 ) 


H(Lk) = LG (L) 

J 

which Is plotted In Figs. 20 and 21 for the von Karraan transverse 
and longitudinal spectra respectively. We shall show later that 
for both of these spectra we have 


H 


- I = -1.666 


(5.43) 


and we see from Figs. 20 and 21 that for values of kL larger than 
about 3 this asymptotic value has, for practical purposes, been 
reached. Thus, for cases where 0 and k L = N LAk ^ SL, we 

may replace Eq. (5.^1) by 


(1) A 


NL^Ak 


H (1) 

[H(5)-H„]d5 + ^ = ^ 


NL^Ak 


+ 


H 

OC 

17 


C5.^4) 


where 


Y 


(1) A 


[H(C)-H_]d?. 


( 5 .^ 5 ) 


(2) 

Let us turn now to I defined by Eq . (5*32) and its large 
T asymptotic form defined by Eq . (5. 33b) — which is 


( 2 ) _ 1 


NAk 


<1 

c 


NT A 

^ £nF(k;L) 


dk 


(5.46) 


where the above limiting form is arrived at by the same arguments 
as were used in obtaining Eq . (5-35). Proceeding along the same 

lines as in the case of we have upon introduction of the 

definitions of Eqs . (5-37) and (5*40), 
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( 2 ) . 


NL^Ak 


N^Ak 


N^Ak 


(Lk)Ldk 


1 


NL^Ak 


N LAk 


N^LAk 




1 


N LAk 


NL^Ak 


N^LAk 


[H^(C)-H^]dC + 


(N -N, )LAk 

oo '■ 2 1 

NL^A k 


1 


NL^Ak 


N LAk 
2 


N^LAk 


[H^(0-H^]d? + 


( 5 .- 47 ) 


where N = N^-N^ was used In going to the last line. Again recog- 
nizing that we are dealing in the integrand with the quantity 
described by Eq . (5*42) ^ we see from Figs. 20 and 21 that for 
lower (dimensionless) wavenumber limits k^L = N^LAk ~ 0 and upper 
limits k L = N LAk larger than about 3L, we may replace Eq. 

(5.47) by ' 


(2) A 1 


NL^Ak 


( 2 ) 

[H"(^)-H^]d? + — = ^ 

NL^Ak 


H' 


( 5 .^ 8 ) 


where 


Y 


(2) A 


[HVC)-H5]dC 


(5.49) 
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From Eqs . (5.29), (5.32), and (5-33) a it follows that the 
large T asymptotic form of the correlation coefficient p can be 
expressed as 


( 1 ) 


P = ~ 




( 2 ) 
i(i) 

CO 


Y 


Cl) 


^«> NLAk 
^ 


NLAk 


H 


(k,-k,)L 




+ 


X 


«= (k^-k^L 


(5.50) 


where the second and third lines are valid approximations when- 
ever we have for the lower wavenumber limit in our likelihood 
equations ( 3 . 32 ) and (3-33) k^ ~ 0, and for the upper wavenumber 
limit we have k ^ 3L-F Notice from Eq . (5.50) and the fact that 
is negative that 


11m 

k2^°° 


p = 1. 


(5.51) 


Hence, when k^ 0 and k^ our estimates of a^L and L have a 

correlation coefficient of unity. This limiting behavior may be 
understood from the fact that if L were known exactly ^ our esti- 
mate of a^L given by Eqs. (3.25) or (3.32) would improve 
indefinitely as k = N 2 Ak increases indefinitely; hence, for 
large enough k^ all of the statistical uncertainty in our esti- 
mate of a^L is a consequence of the uncertalnity in our estimate 
of L — which is reflected in the fact that p -> 1 as k^ 

Let us turn now to evaluation of |A“M given by Eq. (5.30), 
which is required in our expressions for and o\ given by 
Eqs. ( 5 . 27 ) and (5.28). Prom Eqs. (5*30) through (5.33)j we see 
that the large T asymptotic form of |A”^|/N^ can be expressed as 
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1 


(5.52) 


.JCL. 


(a^L) 




Hence, whenever 0 and ^ 3L ^ , we have for practical 

purposes 


I — 1 

II 

1 

< 

M2) _ 

VI 


oo 

\ " / J 


1 j , „ 

L^(a^L)^ I (k^-kj)L 



(kj-k^L 


+ 



1 

(a^L) ^NLAk 




/ ( 2^ 

\y I 

(k^-kpL 


where Eqs . (5-^4) and (5*^8) and the relationship 


NAk = (N2"N^)Ak = k^-k^ 


(I 


have been used In going to the second and third lines In Eq. 
(5-53) • Let us now define 


q3) A 


f>00 

[H(5)-H^]^dC 

0 


.CO 

. 

0 


[HM5)-2H^H(5)+H^]d5 


r 


[H2(5)-Hi,]dC-2H 


[H(C)-H„]d5 


- 2H„Y^^L 




.53) 


.54) 


.55) 
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according to Eqs . (5.^5) and (5.49). Combining Eqs . (5.53) and 
(5.55) with Eq. (5.34), we find that |A“^|/N can be expressed as 


|A'^ I 

N 


1 

(a^L)^L Vt 



(kj-kj)L 


(5.56) 


where y ^ is expressed in terms of y and y by Eq . (5.55) j 
and where Eq. (5.56) is valid for practical purposes whenever 
s; 0 and k^ ^ 3L”^ . 


If we divide the numerator and denominator of Eq. (5.27) by 
N, we can express the large T asymptotic form of as 


o 


2 

1 


( 2 ) 


- 1 


/N 


(5.57) 


where Eqs. (5.32) 
is the variance of 
relative variance 
k^ ~ 0 and k^ L 3L 


and (5 -33b) have been used. Recalling that 
our estimate of we can express the 

of our estimate of a^L for the case where 
using Eqs. (5-48), (5.54), and (5-56) as 


(cr^L) 


H' 


+ 


(3) 


H 


Y 


L 

(IT L 


y(2) ~l 

'(k^-k^)L 

(k^-ki)Ll 

( 2 ) 


1 + 


Y 


/H: 


(k^-k^ )L 


(/D)V/3) 


1 - 


(k^-k^)L 


(5.58) 


Similarly, if we divide the numerator and denominator of 
Eq. ( 5 . 28 ) by N, we can express the large T asymptotic form of 
as 
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1 


(5.59) 


(o^L)^ I A“‘ |/N 


Recalling that is the variance of our estimate of L, we can 
express the velav-ive variance of our estimate of L for the case 
where ~ 0 and k^ ^ 3L“^, using Eq . (5.56) and (5.59)a as 


^ = k 

' (3) 1 


1 

CM 

1 

1 


Y 


1 L 
(3) I 



(k^-k^)L 


(5.60) 


Finally, recognizing that is negative for both the von 
Karman transverse and longitudinal spectral forms, we can express 
the correlation coefficient of our maximum likelihood estimates 
of a^L and L from Eq . (5.50) as 


P = 


1 + 


1 + 


Vh 

* 00 
(k^-k^ )L 


y 


(2) 


/H; 


(k2-k^)L 


( 5 . 61 ) 


which is valid whenever k^ 0 and k^ ^ 3L”^. 

Equations (S.Sd), (5.60), and (5.61) are the main results of 

this section. Each of these three results is the large T 
asymptotic form of its left-hand side that is valid for practical 
purposes whenever h ^ ^ 0 and _> It may be seen from 

Eqs. (3.35), (3.37), (5.37), and (5.40) that for both the von 

Karman transverse and longitudinal spectra , we have 


H 

CJO 


117.97 
70 . 78 


1 . 6667 


5_ 

3 ^ 


(5.6 2) 
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where -6/3 is the exact asymptotto result. The other d-imensi-on- 

(1) (2) (3) 

less constants y j Y ^ ccnd y are defined In terms of the 
function H(^) by Eqs , ( 5 , 4 5 ) ^ ( 5 . 4 9 ) ^ and the first line In 

Eq, (5.55) respectively . The last line In Eq. (5.55) shows that 
only two of these three constants Is Independent . The function 
H(^) = H(Lk) Is defined by Eq. (5.37)^ where (d/dL) InF (k;L) Is 
given by Eqs. (3.35) and (3.37) for the von Karman transverse 
and longitudinal spectra respectively . L Is the Integral scale 
associated with the von Karman spectral forms^ ki and k 2 are the 
lower and upper wavenumber limits used In the likelihood 
equations — see Eqs. (3.32) and (3.33) — and L = VT Is the length 
of the record. 


Coefficient of Variation of Mean-Square Velocity Estimates 

We have shown earlier that, as T our maximum likelihood 

estimates of a^L and L are asymptotically governed by a joint 
normal probability density with mean values equal to the true 
values of a^L and L, and with variances and and correla- 

tion coefficient p, as given above. In particular, Eq . (5-60) 
provides a general expression for the variance of our estimate 
of L (divided by the square of its mean) . 

To get a comparable expression for the variance of our esti- 
mate of where 


(5.63) 


we must consider the ratio of our estimates of a^L and L as indi- 
cated in Eq . (5-63) ■ Since these estimates of a^L and L are 
asymptotically governed by a joint normal probability density, 
the estimate of (obtained by dividing an estimate of a^L 
by an estimate of L) is governed by whatever probability density 
describes the ratio of two (correlated) normal variates each 
having a nonzero mean value. When both mean values are zero the 
resulting density of the ratio is a Cauchy probability density — 
e.g., pp. 153-15^ of Ref. 43 — however, when mean values are not 
zero the resulting density apparently cannot be written in closed 
form — e.g., p. 4ll of Ref. 48. Fortunately though, we may 
obtain an expression for the variance of our estimate of 
using the method of Ref. 49 — even though its distribution cannot 
be expressed in closed form. 
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Merrill (Ref. 49) obtains a series expansion for the vari- 
ance of the ratio of two normally distributed variables which 
converges rapidly (at least in the asymptotic sense) whenever 
the coefficient of variation of the variate in the denominator 
is small — which In our case occurs whenever oj/L^ is small — 
as will generally be the case. Adapting Merrill's notation to 
our case, let 


V 


2 

1 




(5.64 a,b) 


Then Merrill's Iq Is our and Merrill's I is our 
a^. At the top of p. 56 of Ref. 49 ^ Merrill gives 
for the variance of the ratio — l.e., the variance 
mate of — which is 


est-imate of 
an expression 
of our estl- 


Var(a^) = (a^ ) ^ {Vj-2pv^V2+v^+8Vj-l6pvJ v^ 

+ 3vjv2+5p^v^v2+69v^^-138pv^v^ 

+ 15 v^ V 2 + 54 p^ v'J v^+ • • • } , (5.65) 


where we carried through powers of six In the product v^ and V 2 . 

We are primarily Interested in cases where v^ and y\ are 
small In comparison with unity. From Eqs. (5-58) and (5.60), 
we see that v^ and y\ are made small by Increasing the length 
L of the turbulence record. In fact, when Eqs. (5*58), (5.60), 
and ( 5 . 61 ) are substituted into Eq. (5.65) — using Eqs. (5.64 a, 
b) — we see that the right-hand side of Eq . ( 5 . 65 ) takes on the 
appearance of a series in powers of L/1. If we retain only the 
first power in L/L, the relative variance of our estimate of 
becomes 


Var (g^ ) 
(a^)^ 


v^-2pv V +v^+0 
1 '^12 2 



( 5 . 66 ) 


where, according to Eqs. (5.64 a,b), Vj and v^ are given by 
Eqs. ( 5 . 60 ) and (5.58) respectively, and p is given by Eq. 
( 5 . 61 ). The coefficient of variation of our estimate of 

is [Var (g^ )/(a^ ) ^ ]^. 
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The case of most Interest is that for which Ic^ is very 
large. To study this case, we can let (k^-k ^ )L-^oo . For this 
limiting case, we see from Eq. (5 •51) that p = 1. Hence, as 
(k j-k^ )L-><», we have from Eq. (5*56) by retaining only terms 
of first degree in L/L: 




(k “k )L^oo 
V. 2 1 




L a^L 


(3) 




Y 


(k^-k j ( 5 . 67 ) 


where we have used Eqs . (5.64 a,b), (5-58), and (5.60) in going 
to the second and third lines, and where we have retained only 
the first power of L/i in the last line of Eq . (5.67). 


Numerical Results and Discussion 

Let us turn now to numerical evaluation of the above 
described quantities for the von Karman transverse and longi- 
tudinal spectra. We shall treat the longitudinal case first 
since the integrals required for this case are easier to 
evaluate than those for the lateral case. 

The case of most interest is that for which k 2 is very 
large, so we can let (k 2 -ki)L^“ as before. In this case, we 
see from Eqs. (5-58), (5-60), ( 5 . 6 I), and ( 5 . 67 ) that 
( 1 ) 

y is the constant of primary importance — since we have 
already determined that Hoo = -5/3 [Eq. (5.62)] for both the 
von Karman transverse and longitudinal spectra. 

( 3 ) 

The expression we shall use for y ^ is that given by 
the first line in Eq . (5*55): 


Y 


(3) 



0 


[H(0-H^]"d^, 


( 5 . 68 ) 
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where H(C) is defined by Eq. (5.37). For the von Kavman longi- 
tudinal sgeotrum^ we have from Eqs . (5.37) and (3.37): 


H(5) = _ 117.9751 
1+70. 785“^ 


(5.69) 


Hence, H(^)-H^ can be expressed as 


H(^) - H 
^ ^ ' 00 


117.97 
70 .78 


117. 97C^ 
1+70 . 78 ^^ 


which Is of the form 


H(5) - H„ 


a _ a^^ 

^ l+b^2 


a 

b(l+bC^) 


where 


a 

b 


-H = 


00 


5 

3 


(5.70) 


(5.71) 


(5.72) 


and 


b = 70 . 78 . ( 5 . 73 ) 

(3 ) 

Prom Eqs. (5.68) and (5.71) It follows that y can be expressed 
for the longitudinal von Karman spectrum as 


Y 


(3) 




0 


1 

(l+b6^)^ 




(5.74) 


which Is of the form of Eq . (3.251.11) on p. 295 of Ref. 50: 


j (l+bC")"M^ 

0 



(5.75) 
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Hence;, for the von Karman long'itud'inat spectrum ^ we have 


Y 


(3) 



^ 4(70.78)^ 


0.2593. 


(5.76) 


From Eqs . (5.60) and (5*76), we therefore have for the rela- 
tive variance of our estimate of L as (k^ -k ^ )L-»-«> , 


— = 3.856 L/L, (5.77) 


whereas 3 from Eqs. (5.62), (5.67), and (5.76), we have for the 
relative variance of our estimate of as (k 2 ~k^)L^«’, 


^ = 1.714 L/l, ( 5 . 78 ) 

(an" 


where only the first power of h/L has been retained in Eq . 

(5*78). Equations ( 5 . 77 ) and (S. 78 ) apply to the von Karman 
longitudinal spectrum , L is the Integral scale of the turbulence 
and L Is the length of the record. 


Let us now turn to evaluation 
estimates of L and in the case 
spectrum. For this case, we again 
given by Eq. (5.68) where H(^) is 
von Karman transverse spectrum^ we 


of the relative variance of our 
of the von Karman transverse . 

require the evaluation of y 
defined by Eq . (5*37). For the 

have from Eqs. (5-37) and (3-35) 


H(E) = li7.97SUl-l88.75g") 

(1+70 . 785 ") (1+1 88. 75C") 
which is of the form 

H( 5 ) = agUl-cg i ) ^ 

(l+bg")(l+cg") 


(5.79) 


( 5 . 80 ) 
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where 


a = 117.97 

b = 70.78 = I- a 
5 

c = 118.75 = I b . 

From Eqs . (5.80) and (5.81), we can see that 


H 


00 


11m 


H(?) 


ac _ ^ ~ ^ 

be ” b 3 ' 


( 5 . 81 ) 


( 5 . 82 ) 


There fore , 


H(5)-H 

^ 00 


ag^ (1-c^ ^ a 

(l+b5^)(l+c5^) 


a 

b 


(1+bg^) (l+cg^)+b^^(l-cg^) 

(l+bC^)(l+c^2) 


a r l+(2b+c)g^ 

^ d+bC")(l+c^") 


(5.83) 


Combining Eqs. (5.68) and (5.83) gives us the desired expression 
(3) 

for Y in the case of the von Karman transverse spectrum : 


Y 


(3) 


(rV 

.00 

l+( 2 b+c)^^ 

V/ j 


(l+bS^) (l+cCO 



[l+(2b+c)g^]^ 

(l+b5U Vl+o5^)^ 


since the Integrand is an even function of 5. 


(5.84) 
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The integral on the right-hand side of Eq. (5.84) is readily 
evaluated by contour integration using the method outlined on 
pp. 584-587 of Ref. 51. This integration is carried out in 
Appendix I and yields 


7T / 

-Y 

(b+c ) ( 7b ^+c ^ ) 

b (b ^+bc+2c ^ ) 

(c-b)’ ^ 

b/ 

4b'"^ 

J 


(5.85) 


When the exact values for a, b, and c given by the left-hand 
sides of Eqs. (3*41), (3*39), and (3.40) respectively [compare 
Eqs. ( 3 . 39 ) to (3.41) with Eqs. (5.81)] are substituted into 
Eq. ( 5 . 85 ), we find fov the von Kavman transverse spectrum that 

^(3) ^ 0.4437 . ( 5 . 86 ) 

Prom Eqs. ( 5 . 60 ) and (5-86), we therefore have for the relative 
variance of our estimate of L as (k 2 -k^ ) 

^2 

— = 2.254L/L , ( 5 . 87 ) 


whereas, from Eqs. (5.67), (5-82), and (5.86), we have for the 
relative variance of our estimate of as (k^-k^ )L->“'» , 


= 1 . 002 L/L , ( 5 . 88 ) 

(a^)^ 


where we have again retained only the first power in L/l in Eq . 
( 5 . 88 ). Equations (5.87) and (5,88) apply to the von Karman 
transverse spectrum. L is the Integral scale of the turbulence 
and L is the length of the record. 

Equations (5.77) and (5.78), which apply to the von Karman 
longitudinal spectrum, and Eqs. (5.87) and (5.88), which apply to 
the von Karman transverse spectrum, are the main numerical results 
of this section. 

Discussion . It is instructive to compare the values of rela- 
tive variance given by Eqs. (5.78) and (5.88), which apply to the 
von Karman longitudinal and transverse spectra respectively, with 
the values of relative variance of estimates of obtained by 
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squaring and averaging velocity records. The relative variances 
obtained for this latter estimation procedure are 1.73,2 L/L 
and 1.069 L/L [Eqs . (D.28) and (3.21) of Ref. I 8 ] for von Karman 
longitudinal and transverse records respectively. Hence, the 
relative variances given by Eqs. (5-78) and (5 .88) of maximum 
likelihood estimates of — l.e., 1.71^ L/L and 1.002 L/L — 
are only very slightly smaller than the values of 1.732 L/L and 
1.069 L/L for the squaring and averaging estimation procedure. 
Nevertheless, the maximum likelihood method relative variances 
are smaller, as we would expect from the asymptotic efficiency 
normally associated with maximum likelihood estimates. 

The relative variances of 3-856 L/1 and 2.254 L/L for maxi- 
mum likelihood estimates of L for von Karman longitudinal and 
transverse records respectively are, perhaps, of more Interest. 

In contrast to the squaring and averaging procedure used to esti- 
mate a^, reliable estimation procedures and associated variances 
for obtaining the integral scale from velocity records have not, 
in the past, existed. 
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AIRCRAFT RESPONSE EXCEEDANCE RATES 


In Secs. 4.2 and 4.3 of Ref. 19, a series expansion was 
developed for the mean rate of exceedances N^(y) with posi- 
tive slope of a generic aircraft response variable past a 
specified level y. This result was derived for aircraft 
responses to the three component turbulence model of Eqs. 

(1.2) to (1.4) of the present report. However, these results 
are valid only for situations where the three locally sta- 
tionary conditions described in Ref. 19 are satisfied. In 
this section, a new expansion for the mean rate of exceedances 
N+(y) is developed; this new result depends only on the 
validity of the first locally stationary condition described 
by Eq. (1.8). This condition depends on turbulence properties 
only, and is believed to be virtually always satisfied. Thus, 
the results derived herein apply to supersonic aircraft with 
arbitrarily high Mach numbers — as well as to subsonic aircraft 
for which the results of Secs. 4.2 and 4.3 of Ref. 19 should 
apply. However, the derivations contained in the present 
section are considerably more involved than those of Secs. 4.2 
and 4.3 of Ref. 19, and in order to hold the complexity within 
bounds, we have assumed in the present section that the slow 
turbulence component W 0 (t) in Eq. (1.2) is negligible in com- 
parison with the fast component Wj>(t). 


Application of Rice's Formula to Intensity 
Modulated Gaussian Processes 


In the derivation to follow, we shall evaluate Rice's 
expression for the mean number of crossings with positive slope 
per unit time N+(y) of a stationary process past the level y. 

It was shown by Rice [5] on p. 189-193 of the Wax edition 
that, for stationary processes, one has 


N^.(y) 


yp(y,y )dyj 


0 


( 6 . 1 ) 


where p(y,y) is the joint probability density function of the 
aircraft response y and its time derivative y. A derivation 
of Eq. (6.1) also can be found on pp. 45-47 of Crandall and 
Mark [ 25 ]. 
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To evaluate N+(y) from Eq. (6.1), we shall use an exten- 
sion of the methodology developed in See. 2 of this report 
for evaluating the first-order probability density of the 
process y(t). On the first two pages of Sec. 2, we showed 
that the aircraft response y(t) conditioned on the process 
Cf(u) for all -°o<u<_t is strictly Gaussian with zero mean 
value, but generally nonstationary. This conditional re- 
sponse process is denoted by {y (t ) | ap (u) } , -°®<u£t. Therefore, 
the joint probability density function of y(t) and y(t) con- 
ditioned on cf(u) for -“<u£t is a joint normal density 
function with zero mean values for y(t) and y(t). Let us 
denote this joint conditional density function by pCyjylcp), 
where Of denotes the infinite dimensional "vector" ap(u)”^ 
-“<u£t, as described on the third page in Sec. 2 of this 
report. We therefore may formally express the unconditional 
joint probability density of the aircraft response and its 
first derivative as 


p(y,y ) 


.oo 

p(y 5 y1 5 

0 


( 6 . 2 ) 


where p(up) denotes the probability density of the infinite 
dimensional vector qp = Cp(u), -“<u£t, and where the symbolic 
integration is taken over this same infinite dimensional space 
as described on the third and fourth pages of Sec. 2 of this 
report. Substitution of Eq. (6.2) into Eq. (6.1) yields 


N^(y ) 


-00 -00 

y I p(y,y I qp)p(Sp)da^dy 

0 


(6.3) 


When Up =: ap(u), -°°<u<_t is specified, the joint density 
of y and y is a joint normal density with zero mean. values in 
y and y [Ref. 39, pp. 1^7, 1^8] - l.e., 


p(y.ykp) = p(y,y| cr^,a|,ii^^) 


(cr?y^-2]j .yy+a^y^) 

y yy y 

2 (a ? -|i ^ . ) 

y y yy 


2iT(a^a?-y ^ . ) 

y y yy 


(6.4) 


138 




each of which is a stochastic function of time t that is deter- 
mined by the behavior of the stochastic modulating process 
Ql'Cu) over the time interval -°°<u£t. Thus, the expectation 
operations in Eqs . (6.5a) to (6.5c) are taken over variations 
in the (stationary Gaussian) modulated process (z(t)} in our 
turbulence model of Eqs. (1.2) to (1.^). [Recall that, in 
this section, we are assuming Wg(t) =0.] We have indicated 
that these expectations are with respect to the process {z(t)} 
by placing the subscript "z" after the expectation operator 
"E" in Eqs. (6.5a) to (6.5c). 

Noting that Eq . (6.2) expresses the mathematical expec- 

tation of the response joint conditional density p(y,y|c^) 
with respect to the modulating process (af(u)}, -°°<u<_t, we 
can also express Eq . (6.2) as the expectation 

p(y,y)=E{p(y,yla)}, (6.6) 

f 

which when combined with Eq. (6.1) yields 


N^(y) 



0 


yE^ {p(y,y I C^)} dy 


(6.7) 


which we shall now proceed to evaluate in terms of measurable 
metrics of the modulating process {a^(t)l. 
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Series Expansion of Conditional Joint Probability 
Density of Aircraft Response Displacement and Its Derivative 

As in Sec. 2, we are Interested here in excitation 
processes Wf(t) = Of(t) z(t) where typical fluctuations in 
a£.(t) are not more than about one-third* of the mean value 
of Of(t). Different sample functions ap(u), -<»<u£t give 
rise to different values of the three parameters ay(t) 

and ]jy^(t) in the probability density, Eq . (6.4), 

as may be seen from Eqs. (6.5a) to (6.5c). Hence, in a 
manner analogous to the approach used in Sec. 2, we require 
here a truncated Taylor’s series representation of the joint 
normal density, Eq. (6.4), in the three-dimensional para- 
meter space a?., and *^he expansion will be centered 

about the point defined by the mean values a?, and u . 

^ • y y ^yy 

of the variables a^(t), a^(t), and p .(t): 

y y yy 


*The approximate upper limit of one-third for the typical 
fluctuation 5a^(t) relative to the mean Of is arrived at as 
follows. Denoting expected values by overbars and fluctua- 
tions by delta, we have 


a = G + 5a 


(a) 


hence 

G^ = (g)^ + 2g5g+(5g)^ , (b) 

and 

a ^ = (g ) ^ + (6g ) 2 , ( c ) 

since 5g = 0 by definition. Therefore, from (b) and (c), we 
have 

G^ - G^ = 2g6g + (6g)^ - (6g)^ 

2g6g , (d) 

since, by assumption, we have 5g << g. From (d) it follows 
t hat 
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(6.8a) 


- E {a^(t)> 
y y 

a? = E {a^.(t)} 

y y 


and 


yy 


= E 


{y • ( t ) } 
^yy 


( 6 . 8b ) 


(6.8c) 


where we have again used the subscript on the expectation 
operator to denote the process that the expectation is taken 
with respect to. 


E{(a2-a2)2} ^ i|(a)2E{(5a)"} . (e) 

From (c), we also have when 6a << a 

^ :: (a)^; (f) 

hence, combining (e) and (f) yields 

E{(g^^)^} ^ ^ E{(6g)M ^ 

(a^)^ (a)^ 

where E{(6a)^} is the variance of 6a since 6a e 0. Applying 
(g) to ay, and using the definition of y given by Eq. (2.27a), 
we have 


'E{(5ay)')' 


rE{(CT^-a0^11 
y y 

h 

^ 1 

(~ )" 

L ^ y J 


L ' y J 

/4y 


(h) 


Moreover, from Figs. 9 to l6 , we see that the largest value 
of Y that gives "^ood" accuracy is y = 2. For this value, we 
have ( ijy )-'! E 8 2 = 0.35^ which is slightly more than one- 
third. Therefore, since typical relative fluctuations in 
ay(t) should never be larger than those of a^(t), we conclude 
that whenever typical fluctuations in a^(t) are not larger 
than one-third of the mean value of ap(t), our methods should 
provide reasonably accurate results. 


I4l 



In order to evaluate a common form of the multi- 
dimensional Taylor’s series [e.g., p. 338 of Ref. 35], we 
require an expression for the square of a trinomial: 

(x+y+z)^ = (x+y)^ + 2(x+y)z+z^ 

= x ^+y ^+z ^+2xy+2xz+2yz . (6.9) 


By combining Eq. (6.9) with Eq. (57) on p. 338 of Ref. 35, 
we obtain a Taylor’s series representation of the right-hand 
side of Eq . (6.4) about the point ' Taking the 

the expected value of this representation with respect to 
the process ap(t) yields in straightforward fashion 


E^ {p(y,ylgf)} = p(y,yla^, a|, y^^) 


+ E {(o"-aP}p^^»°’°ky,y|o^ a5. y .) 


f 


y y 


y y yy 


+ E„ {(cri-ap}p(°>l’°ky,y|a% a?, y .) 


y y 


y y yy 


+ E {(jj .-y •)}p^*^’^’^^(y 3 y|o'^, y .) 

yy yy wjji yj yj ^yy/ 

1 ^ r / 2 _rr ^ 2 1 ^ ( 2 , 0 , 0 ) / ^ I 7TT 




y^ y. ^-yy 


+ ^ E {(p^-a?)^}p^'^’^^^^(y,yla^, o^., y,, •) 

^ O’ 'tT Tt' y \ y \7 ' trTT 


y y 


y yy 


+ i E {(y .-TrT)"}p^°’°’^^ 
2 ' "^yy ^yy ^ 


(y.yla^, o^. , y .) 

V J i J 1 y 5 ^-yy 


f 


y ’y' "y ’y" ^ ' y' ' y' ^ yy 


+ E„ { (a5-a,^J (y ,y I a% y .) 


+ E { (a ^-a ^ ) (y . -y . ) } p 

y y '^yy yy 


(1,0,1) 


(y.yla^, a^, y .) 

V j , j I y 5 y > Hyy 


+ E„ {(a?-a|)(y„.-y .)}p‘-°>^’^4y,y|a,% oL 


-f y y yy yy 

+ higher order terms. 


y" yy 


( 6 . 10 ) 
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where Eq. (6.9) 
In Eq. (6.10). 


was used in writing out the last six terms 
In Eq. (6.10), we have used the definition 


p(ijJ sk) 




y 


y • ) 
yy 


a(a^)'"8(a?)^9(y .)^ y y yy 


y 


yy 


2 2 
a = a ^ 


2 2 
a f = a 7 


yy 


= y 


yy 

( 6 . 11 ) 


is the joint normal density given by 


where p(y,y|a^, , y •) 

y y yy 

the right-hand side of Eq. (6.4), and where the partial 
derivatives are evaluated at the expansion center defined 
by Eqs . (6.8a) to (6.8c). The partial derivatives evaluated 

at this expansion center are not random variables. Thus, 
the expectation operation with respect to the process Cf(t), 

E {p(y 3 yla„)}, in the left-hand side of Eq . (6.10) yields 

Q f "“I 

the expectations {•••} of the various expansion coefficients 

in the right-hand side of Eq . (6.10). In Eq. (6.10), we have 

Included all terms with partial derivatives of order 1+J+k = 2 
or less, which contain all terms with powers of two or less 


in products 


( P \7-\*7-“y yy ) 


of the random variables (Oy-Oy), (a^-a|-), 

. y-^.yy/ . It Is an implicit property of the Taylor’s series 
expansion that the terms written out in the right-hand side 
of Eq. (6.10) are the most significant terms ~ provided that 
fluctuations in a^(t) about its mean values are not too large 
as indicated earlier. Since, we have 


and 


E {(a"-a")> = 0 

y y 


aa|-a|)} = 0 


{(w .-U .)} 
. yy yy 


= 0 


( 6 . 12 ) 
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the first three "correction" terms to p(y,y|aya a?, 

in the right-hand side of Eq. (6.10) are identically zero. 
We shall now proceed to develop expressions for the six 
remaining correction terms written out in Eq. (6.10). 


Evatuati-on of devi-vatives of foint Gaussian prohahitity 
density. The somewhat tedious job of evaluating the six 


partial 

required 
report . 
values 


derivatives p*" ^ (y sy 1 

in Eq . (6.10) is carried out in Appendix J 
The resulting expressions evaluated at the 


a! , 
y 




l+j+k = 2 


of this 
expected 


of a?, and ]j . for use in Eq. (6.10) are given 

y y yy 


by Eqs . (6.15) to (6.20). 
we have used the fact that 


In evaluating these expressions 


yy 


EE {y .} = 0 

yy 


(6.13) 


which is shown later by Eq. (6.7^). Equation (6.l4) gives 
the joint density of y and y also evaluated at the expected 
values of and y .: 

y y yy 


p(y =y I Oy 

p(2,0,0) 

p(0,2,0) 

p(0,0,2) 


(y^yky 

(y^yky 
(y ,y I cTy 






2 a 


y 


2 a? 


2tt (a ^ a ^ ) 

y y 


2 > ^ 






2 4 

y y 


3 - 6 ^ + 

(cr^)^ 

y y 




Map- 


3 - 6 + 


(7^ (a^^ 

y y 


y • ) 
yy 


p 


2 2 
a a . 

y y 




1 - 3 = 1 - 


iL 


y 


a? 


(6.14) 

(6.15) 

( 6 . 16 ) 

(6.17) 
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(6.18) 


(l^ljO)^ .1 2 ^2 

y" y yy 
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i|a' c? 

y y 


(I.O5I)/ .1 2 2 \ 

(y 5y 5 cr: , y • ) 
' y y yy 
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- 3 
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(at) 
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(6.19) 


P^°’^’^^y,y|cf^, a;, Uy^) = 


y 


2(a^)'"^(a!)'/^ (a!)"'^ 


y y 


y 


y 




- 3 


y 


1. 

2 \'^ 


iap 


(6.20) 


In each of the above six relationships, p is the joint normal 
density function given by Eq . (6.l4); furthermore, in arriving 
at each of these six relationships, we have used the fact that 
y . = 0 as indicated by Eq . (6.13)- 


Series Expansion for Aircraft Response Exceedance Rates 

According to Eqs . (6.7) and (6.10), we must now multiply 

each of the above seven expressions, Eqs. (6.14) to (6.20), 
by y and Integrate each resulting expression with respect to 
y from zero to infinity. Examining Eqs. (6.14) to (6.20), we 
see that these integrations require evaluation in the follow- 
ing expression for values of n from 1 to 5: 



0 


( 6 . 21 ) 


These integrals are easily evaluated with the aid of Eqs. 
(3.461-2) and (3.461-3) on pp. 337 of Ref. 50 yielding: 
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I(n) 



and 


I(n) 



n = 2,4,6 


3 


where in Eq. (6.23) 


we have used the definition 


(n-1) ! ! = !• 3*5* • • (n-1) , 


( 6 . 22 ) 


(6.23) 


(6.24) 


and in Eq. (6.22) we have used the usual definition of a 
factorial. Evaluation of Eqs. (6.23) and (6.24) for n from 
1 to 5 gives 


1 ( 1 ) 



1(2) = 1(1) 

1(3) = 2 a| 1(1) 

1(4) = 3/|y'^(aT)''^=‘l(l) 

and 

1(5) = 8 (^)U(l) . 


(6.25) 

(6.26) 

(6.27) 

( 6 . 28 ) 

(6.29) 
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When Eq . ( 6 . 10 ) is substituted into Eq. ( 6 . 7 ), we see 
that we must multiply each of Eqs . ( 6 . 14 ) through ( 6 . 20 ) 
by y and Integrate the resulting expressions from zero to 
infinity. When these integrations are carried out with the 
aid of Eqs. ( 6 . 25 ) through ( 6 . 29 ), we find that the mean 
threshold crossing rate with positive slope N_j_(y) can be 
expressed as 


N^(y) - 



_ JL 


20 . 






3 _ 6 ^ 

a.l (a!)" 
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1 ^f 


{ (a ) ^} 

y y 
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{ (y . ~y . ) ^} 

. yy yy 


2 2 
O G . 

y y 


1 - 
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y y y y 


1 - 


2 2 
a ^ a . 

y y 
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-(i)' 


{ (a^-a^ ) (y . -y . )} 

. y y ^yy *^yy 

«/ tX 


2 ^ ^ 






+ higher order terms. 


( 6 . 30 ) 
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In arriving at Eq. (6.30) , we have made use of Eq. (6.12) 
and the fact that 


.00 


0 


y p 


( 0 , 1 , 1 ) 


(y,y 




y . ) dy = 0 . 

yy 


(6.31) 


The "higher order terms" In Eq. (6.30) arise from the terms 
with the same label In Eq. (6.10). As mentioned earlier, 
the "slow" turbulence component Wg(t) In the turbulence model 
of Eq. (1.2) has been assumed to be zero In arriving at Eq . 
(6.30)* We also have used the fact that tyy ^ 0 In arriving 
at Eq. (6.30). Later on in this section^ we shall show that 
the last displayed term in Lq, ( 6 . ZO ) is identially zero^ 

[Eq, (6.71)[\. 

The problem of finding expressions for the expansion 
coefficients In Eq. (6.30) will be addressed next. 


Expression for Exceedance Rate Expansion Coefficients 

Here, we shall derive a general expression for the 
expansion coefficients E^ {•••} In Eq. (6.30). We shall 

then show how the various Individual expansion coefficients 
In Eq. (6.30) can be computed for an arbitrary aircraft 
modeled as a linear time-invariant system. 

All three conditional response variables defined by 
Eqs . (6.5a) to (6.5c) can be expressed by the single quantity 

yjk = ” E^{yf t)y^(t) laj,(u)} , - <”<u<t, 

j + k = 2, (6.32) 

where superscripts j and k denote powers of y(t) and y(t) 
respectively, and Eg { • * * | cr^( u ) } , -“<u£t denotes expectation 
with respect to the process {z(t)} conditioned on the process 
(ap(u)} for -«><u£t, where z(t) and Cf(t) are components In 
our turbulence model of Eqs. (1.2) to (1.^). Thus, by 
comparing Eqs. (6.5a) through (6.5c) with Eq . (6.32), we 

see that 


1^8 



2 — 

a = y ., 
y Jk 

for 

J = 2, 

k = 0 

(6.33a) 

a? = y 
y Jk 

for 

j = 0, 

k = 2 

(6.33b) 

y • = y -1 
yy jk 

for 

J = 1, 

k 1 

(6.33c) 


Using the notation of Eqs. (6.32) and (6.33 ) 5 we can now 
express any of the expansion coefficients In Eq. ( 6 . 30) In 
the general form 


E„ Um. 


-y 


J’k’ 


t ) (y 


"k 


!t~y ^ ^ 


= (t)y^ (t)|c^] - E^Cy*^' (t)y^ (t)|q^]j 

X (E^[y^"(t)y‘""(t)|a^] - EjyJ"(t)y‘""(t)|g^] jl (6.3^) 

where we have used to denote conditioning on a£>(u), -<»<u<_t 
as before, and where^the overbar denotes an expectation with 
respect to the process (af(t)}. The Individual expansion 
coefficients in Eq . ( 6 . 30) that are particular cases of 

Eq. (6.3^) are obtained using the values of j', k’, j", and 

k” shown in Table 1. The last entry, E^ {(a^-a?)(y -y . ) , 

’ yf y y yy yy ’ 

is not required for use in Eq. (6.30) because its "multiplier” 
was shown to be identically zero. 

To obtain a general expression for yji<;(t) defined by 
Eq. (6.32) that covers all three cases listed in Eq . (6.33) , 
let us now define the generalized "instantaneous cross- 
correlation function" of the aircraft impulse response as 


J k 


h . (t - i) 

J ^ 


hk(t + ^) 


(6.35) 


where hn*(t) and hj^(t) are tabulated in Table 2 for the three 
cases listed in Eq. (6.33). Functions h(t) in Table 2 are 
the time-derivatives of the aircraft impulse response function 
h(t). Since convolution of h(t) with the input yields y(t) 
as indicated by Eq. (1.21), it follows directly from Eqs. 
(K.l4) and (K.15) in Appendix K that Table 2 contains the 

appropriate definitions of h.(t) and h, (t). 

J 
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TABLE 1. VALUES OP THE PARAMETERS IN EQ. (6.34) THAT YIELD 

THE INDIVIDUAL EXPANSION COEFFICIENTS IN EQ. (6.30) 


















To apply the results of Appendix K to the computation 
of y., (t), we further let 


Xj.(t) = x^(t) = w^(t)|a^(u) = a^(t)z(t) 


(6.36) 


- «><U<o° 


in Appendix K, where Wf(t), a^Ct), and z(t) are components 
of our turbulence model of Eqs . (1.2) to (1.4), and where 
the conditioning operation in Eq . (6.36) indicates that crf(t) 
is to be treated as a known function. With the interpreta- 
tions provided by Eqs. (6.35) and (6.36), the instantaneous 
cross-spectrum input-response relation of Eq . (K.21) becomes 


*y.y, 


2f) 


(f ,t-u| g„)$, . (f,u)du 

Wf 1 


(6.37) 


where 5> (f,t|a„) is the instantaneous cross-spectral 

density of Pj]^(t) [defined by Eq. (6.32)], (f,t) is the 

3 k 

Fourier transform with respect to t of Eq . (6,35) as inter- 

preted by Table 2, and (f,t|g„) is the Fourier transform 
with respect to t of ^f 


(T,t 


w 


= 


a^(t 


- 


a^(t 






( 6 . 38 ) 


where (j) (t) Is the autocorrelation function of the component 
z(t) or the turbulence model of Eqs. (1.2) to (1.4), and 
Cf(t) is to be considered as a known function in Eq . (6.38). 

When we apply the locally stationary approximation pro- 
vided by the first term on the right-hand side of Eq. (2.45) j 
we have 



( f ,t I c^) 


c^(t )$^(f ) , 


(6.39) 


where ^z(f) is the Fourier transform of (}> (t). By substi- 
tuting Eq. (6.39) into Eq. (6.37) and integrating over all 
f, we obtain 
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_0O ^ ^ 

f oo ,00 

a^(t-u) j h 5 

L J k 


(6.40) 


where we have Interchanged orders of integration in going to 
the second line, and where we have used the fact that 
0 (f,t|g^) is the conditional Instantaneous cross-spectral 

density of {y^(t)} and (y (t)}, j + k = 2, whose integral 
over -oo<f<co yields Pjv(t) — see Eq . (6.32) and the general 
property of Eq. (K.ll). By extension of Eq . (2.54), let 
us now define 


Y 


h .h. 

J k 







where as noted above, we have 


\.h (fv) = 
J k 


h^.(t - |) e (6.42) 


where h.(t) and h^(t) have the int erpret at ions listed in 

Table 2'3 Substituting Eq . (6.4l) into Eq . (6.40) yields the 

desired form for y (t): 

Jk 


Pjj^(t) = cr^(t-u)y 


h .h, , z 
3 k’ 


( u ) du . 


(6.43) 


In particular, taking the expected value of y., (t) with respect 
to the process a^(t) gives ^ 
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^h.h, 5 

J 


hence, from Eqs . (6.43) and (6.44), we have 


- ^jk 


C°f(t-u)-Of]Yh.h, ,z('")du 

J k’ 


We recall now that each coefficient In the left-hand 

column of Table 1 Is of the form of Eq. (6.34); thus, each 

coefficient Is the expected value with respect to {Qf(t)} 

of cross-products of the form of Eq. (6.45). Furthermore, 

we see that Eq. i_6.45) Is the convolution of the stochastic 

function [af(t)-a|>] with the "system charact erlzat Ion" 

Y, , (t). Thus, we can use the result, Eq . (K.26) 

' h . h, , z ^ ^ 

of Appendix K, appropri ately Interpr eted , to obtain a general 
formula for (y ^ ^ , -y ^ , ) ( t j ^ ) } . To apply Eq . 

(K.26) to this situation, we require the Fourier transform 

of y, , ( t ) : 

' h .h, , z ^ 

J k" 


, . A 

^h .h , z 
J k’ 


/ , ^ 12TTVt 

^h.h, dt 

J k’ 


Thus, applying Eq . (K.26) to the present situation gives 

’k'“^j 'k’ ^ "k"^^ 


* 2(v)y^ h z*^“^^^h h ’ 


where $ z(v) Is the power spectral density of the process 


1 



{o^(t)-a^} which has no d-c value, the functions Yj^ (v) 

are defined by Eq . (6.46) with primes added as appropriate, 
and the superscript asterisk denotes the complex . conj ugate . 
The negative arguments in Eq. (6.47) arise from the fact 
that we have used a positive exponent in Eq. (6.46), whereas 
the transforms Hj(f) and H^(f) in Eq. (K.26) are defined as 
in Eq. (1.9) with negat ive^exponent s . Equation (6.47) is 
the main result of this subsection. 


Expressions For System Characterizations 

Here, we examine the system characterizations with respect 

to the process {z(t)} that are defined by Eqs. (6.4l) and 

(6.46). First, we note that y, ^ (t) also can be expressed 

h • n , z 

as J k 


Y 


h . h 
J k 


(t ) 


-CO 

$ (f) Re , (f,t)df 

z h .h ' ’ " 


(6.48) 


where Re (f,t) denotes the real part of ^ (f,t). To 

k j k 

show this, we note from Eq . (6.42) that ^ (f,t) can be 
expressed as J k 


J k 


h^. (t - ^)kj^(t + ^) cos ( 27rfT )dT 


-1 


h . ( t 
J 


2)h^(t 


+ ^) sln( 27rfT ) dx 


(6.49) 


Since hj(t) and h]^(t) are real, the cosine Integral in Eq . 
(6.49) IS real, whereas i times the sine integral is 
Imaginary. Moreover, the sine integral is an odd function 
of f. Since even function of f, the contribution 

of the sine integral in Eq. (6.49) to the integral ^ 



in Eq. (6.4l) is zero. We are left with the contribution 
of the cosine integral, which is the real part of (f,t), 

as indicated by Eq. (6.^8). Since real, it follows 

that Yn (t) also is a veal function of t* 

J k’ 

From the fact that y, , (t) is real, it follows from 

Eq. (6.i)6) that 


'i'h.h. 


f- 


(6.50) 


where the asterisk denotes the complex conjugate. Hence, 
Instead of Eq . (6.47), we also can write 


{(li 


j ’k’ ^j ’k 


)(y 


j"k" ^j’'k' 


)} = 


f 


(v)y 


h 


t^ki 


(v)y£ 

z ^ ^ 'h 




(v)dv 


(6.51) 


Moreover, since both quantities Hii^(t)-iij^ in Eq . (6.51) are 
real as can be seen from Eq . (6.45) and the fact that 

Yv, b r^(^) is real, by taking the complex conjugate of both 
k 

sides of Eq . (6.51) it follows that we also have the second 

alternative form 


'k' 'k' ^ "k" '^J"k"^^ 


$^2(v)yS 


(\')Yu >, „(v)dv. 


f 'hj„h^„,z 


(6.52) 


where we have used the fact that the power spectral density 
^>^ 2 (v) of {c^(t)-ap is real. 
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We see — e.g., from Eq . 
system characterization ^ 

j 

is the quantity required_for 

coefficients E {(y.,, ,-y.,, 

j ’k' j ’k 


( 6 . 51 ) — that the transformed 
(v) defined by Eq . (6.46) 

computation of the expansion 
, ) (Pj * Inserting 


Eq. (6.4l) into Eq. (6.46) and interchanging orders of 
integration yields 


^h.h, = f 

1 k 

^ —00 


, J 


( 6 . 53 ) 


- |-)df 


(6.54) 


where we have used Eq . (K.23), and Eq. (K.8) applied to h.h, . 
Functions Hj(f) and H].^(f) in Eq. (6.54) are the Fourier ^ 
transforms of hj(t) and h^(t) as defined by Eq. (1.9) where 
individual interpretations of h^(t) and hj^(t) are given in 
Table 2. 


According to Table 2 

Ti (^^) which are y, , 

' h .h, , z ^ 'hh , 

J k’ " 


5 

z 


we 

(v) 


require three forms 


Y 


hh 3 


( V ) , and y 


hh 5 z 


of 

(v). 


We can immediately write out the first of these three forms 
using Eq . (6.54) as 




.00 

|-)H*(f+ |)df. 


( 6 . 55 ) 


which we have discussed earlier — see Eq. (2.62). 

In order to put ”^hh z^^^ form of Eq . 

( 6 . 54 ), we require the Fourier transform of the time derivative 
of h(t). Let us denote this Fourier transform by placing a 
dot over H(f) — i.e., we def-ine H(f) as 


H(f) 



h ( t ) e 


-IfTTft 


dt 


— 00 


(6.56) 
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The Fourier mate to Eq . ( 6 . 56 ) is 


h(t ) 


j-OO 

^00 


H(f . 


( 6 . 57 ) 


Furthermore, the Fourier mate to Eq . (1.9) is 


( 00 
— 00 


( 6 . 58 ) 


Differentiating both sides of Eq . ( 6 . 58 ) with respect to t 

gives 


h (t ) 



— 00 


12TrfH(f ; 


(6.59) 


hence, comparing Eqs . (6.57) and (6.59) yields 

H(f) = 127TfH(f) , ( 6 . 60 ) 

which relates the Fourier transform H(f) of h(t) to the 
Fourier transform H(f) of h(t). 

Using Eqs. (6.5^) and ( 6 . 6 O), we now can express 
Y,-,* (v) as 

$^(f) (-l)2ir(f+ |-)H*(f+ |■)(l)2^T(f- |)H(f- |)df 


nn , z 


7 . . = 

hh’ 


= 4 7T ^ 


f 

f^<D^(f)H(f- ^)H»(f+ ^)df 


- TT^V^ $ (f)H(f- 
z 

—00 


|)H*(f+ |)df . 


( 6 . 61 ) 
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From Eq. (6.55) a we see that the Integral in the second term 
Is hence, expressed as 




= 4 


7T 


f 


0^(f)H(f- |)H*(f+ ^)df-Tr2 


2 ~ 
V Y 


hh,z 


(v) . (6.62) 


However, the autocorrelation function of the process (z(t)} 
whose sample functions are the time derivatives of z(t) is 
the negative of the second derivative of the autocorrelation 
function of (z(t)} — l.e., -c{)”(t) [see, for example, p. 21 of 
Ref. 37]. Furthermore, by differentiating twice with respect 
to T the relationship 


<1>^(t) = 


^>^(f)e^^'^^'^df 


(6.63) 


one finds 


(P'Ut) = 


47r^f^$ (f)e^^'^^'^ 


df , 


(6.64) 


from which it follows that 4tt f^^> 2 (f) is the power spectral 
density of the process {z(t)}. We therefore can express 




( 6 . 65 ) 


where we have defined 


.CO 

= 4tiH |-)H*(f+ |-)df 


( 6 . 66 ) 


Finally, we consider Yj^j^ combining Eqs. (6.54) 

and (6.60), we have in this case 
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= f 5-)(i)27r(f- |-)H(f- |-)df 


- 12ty 


f$^(f)H(f- ^)H*(f+ ^)df 


- ITTV 


$^(f)H(f- |)H*(f+ |)df 


(6.67) 


However, from Eqs. (35a) arri (35b) of Ref. 3^, we see that 

[H(f- j)H*(f+ ^) ] is necessarily an even function of f since 

$]^(fjt) is an even function of f. Furthermore, <l> 2 ,(f‘) also 
is an even function of f. Hence, we always have 

.00 

f<l^(f)H(f- ^)H*(f+ ^)df = 0 ( 6 . 68 ) 

^OO 

since the Integrand in Eq . ( 6 . 68 ) is an odd function of f. 
Thus, combining Eq . ( 6 . 67 ) with Eqs. (6.55) and ( 6 . 68 ), we 
have 




-iTTVy, . (v) 

hh, z 


(6.69) 


which is the desired expression for (v) . 

^ ' hh , z 


Summary. Equations ( 6 . 55)5 ( 6 . 6 5 ) - ( 6 . 66 ) , and (6.69) 
provide expressions for Yhh,z^^)’ ^hh "^hh 

a function of the power spectral density of the process 

{z(t)} of the turbulence model of Eqs. (1.2) through (1.^) 

and the complex frequency response function H(f), Eq. (1.9)> 

of the relevant aircraft response variable. These expressions 

for y,, (v)), (v), and y ^ (v) are to be combined with 

i^jn ^ z in in ^ 2 inm y 2 

Eq. ( 6 . 51 ) ;to obtain the__various coefficients 

Ea ^ 'k' ^ applying the rules in 
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Tables 1 and 2. By carrying out this procedure, we shall 
now show that the second to the last coefficient in the 
left-hand column of Table 1 is identically zero. 

According to Tables 1 and 2, we have from Eq. (6.51) 
for the second to the last coefficient in the left-hand 
column of Table 1: 




y y 


yy yy 


= ITT 


(V ) 


^hh, z ^ 


dv, (6.70) 


where we have used Eq. (6.69) in going to the second line. 

The integral in the second line of the right-hand side of 

Eq. (6.70) is real; hence, the entire right-hand side is 

imaginary. However, the left-hand side of Eq. (6. 70) is 

real. It follows that the integral in the second line of 

the right-hand side must be zero. This property also follows 

from the fact that the integrand in the second line of Eq . 

(6.70) is an odd function of v since $ z(v) and | y, . „(v)|^ 

cT ^ nn , z 

are both even functions of v. The even property of 
I Yj^h follows from Wiener's theorem and the fact 

that y, , ^(t) is real. Equation (135) of Ref. 3^ is a state- 

n ii ^ ^ 

ment of Wiener’s theorem, where we remind the reader that 
P^(t) is an even function of t. We therefore have 


{(af-a!)(y„,-y„,)} e 0 


y 


yy yy 


(6.71) 


The validity of Eq . (6.7I) also could be argued on physical 
grounds from the fact that it is the coefficient in Eq. (6.30) 
of an odd function of the response variable y. 

Furthermore, we shall now show that 

y.=E{y.(t)}=0 ( 6 . 72 ) 

yy yy 
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as was stated earlier without proof by Eq. (6.13)* Prom 
Eqs. (6.^^) and (6.46), we have 


■jk 


= ol. 


Y 


h.h. ,z 


( t ) dt = a 


"^h.h ,z 
J k" 


( 0 ) 


(6.73) 


Hence, for j - 1 and k = 1 , we have according to Table 2 and 
Eq. (6.73), 


77 - ^ 2 


yy 


= a 


f ^hh,z^°^ ° 


(6.74) 


where the right-hand equality follows from Eq . (6.69). 


Final Expression For Aircraft Response Exceedance Rates 

Incorporat ing the result, Eq. (6.71) into Eq . (6.30), we 

obtain now our final expression for the mean rate of exceed- 
ances N_j_(y) of our aircraft response variable past the level 

y * 


N^(y) - ^ 



2a 


1 - 


E^ {(a?-a?) } 

y y 


8(an 
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0 -^ " y y 
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3 - 6 L + 
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yy yy 
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2 ^2 
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+ higher order terms. 


(6.75) 
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where the "higher order terms" are the same as those in 
Eq. ( 6 . 30)3 which arise from the higher order terms in 
Eq. (6.10). Neglect of these higher-order terms should 
not lead to substantial error provided that typical fluc- 
tuations in the modulating process Cf(t) are not larger 
than about one-third of the mean value of cf(t). The 
turbulence model for which the result Eq . (6.75) applies 
is that of Eqs. (1.2) to (1.4) 3 where here the "slow" com- 
ponent Wg(t) is taken to be negligible in comparison with 
the "fast" component w^(t) = cj^(t)z(t). 

To evaluate the parameters and coefficients in Eq. 
( 6 . 75)5 we require 





_co 




and 


a? = f $ (f)|H(f)| Mf 

y J w ' ‘ 

—00 ^ 


(6.76) 


(6.77a) 



—00 


f^O (f)|H(f) 


df 


(6.77b) 


where $ (f) is the power spectral density of the turbulence 

w^ 

process {wp(t)}3 H(f) is the complex frequency response of 
the aircraft response variable of Interest as defined by 
Eq. ( 1 . 9)5 fi(f) is the Fourier transform of the time 

derivative of the Impulse response function of the aircraft 
response variable of Interest. Equation (6.77a) is a direct 
consequence of Eq . ( 1 . 21)3 and Eq. (6.77b) follows from Eq. 
( 6 . 60 ). Furthermore 3 from Eq . (6.51) and Tables 1 and 2, 
we find that the coefficients in Eq . (6.75) can be expressed 
as 
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( 6 . 78 ) 






", '<»|-«s>‘> ■ ' 




(6.79) 


E { (y . -y . ) ^} = 
'^yy '^yy 




( 6 . 80 ) 


and 


E {(a"-a")(a2-a^} 

y y y y 



— oo 


^^ 2 {v)y ^(v)y* (v)dv 

hli,z 


( 6 . 81 ) 


$ z(v) is the power spectral density of the process 


{af(t)-a|‘}. Two methods for computing this spectrum are 
given on pp . 79-83 of Ref. 19. System characterizations 
(with respect to the process {z(t)}) y.. (v), y^^ (v), 

and given by Eqs . (6.55), (6 . 65 )-( 6 . 66 ) , and 

( 6 . 69 ) respectively. ^z(f*) the power spectral density 
of the process {z(t)} which can be taken as the normalized 


-CO 


spectrum $ (f)/ 

w n 


$ (f)df because of the locally stationary 
w ^ 


assumption Eq. (1.8a), and the normalization assumption, 
Eq. (1.4). 
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Limiting cases of the result given by Eq. (6.75) can 
be studied in a manner analogous to that used in Secs. 2.5 
and 2.6. It can be shown that the result given by Eq. 
(6.75) reduces to the appropriate result in Sec. 4.3 of 
Ref. 19 when fluctuations of a^(t) are assumed to be 
negligible over the duration of the aircraft impulse 
response h(t ) . 
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APPENDIX A 


DERIVATION OF CORRECTION TERMS TO THE GAUSSIAN PROBABILITY 
DENSITY FUNCTION FOR USE IN EQUATION (2.12) 

Equation (2.9) expresses the correction terms in Eq. 
(2.12) as derivatives of the Gaussian density function: 




,k 


2 a 




y 


!7ra 


y 


y 


2 2 
0 =0 

y y 


(A.l) 


where 


(y I ^ 


2a ^ 

e y 


’ y 

Differentiating Eq. (A.l) with respect to Oy yields 


(A. 2) 


|o^) = ^ 


/2tt 


1 3/2 2a ^ 


y 


/ 2 \ y^ r 2\'-2 

+ (ap e y ^ (ap 


p(y kf) 

2 a^ 

y 


- 1 


y 


(A. 3) 


in agreement with Eq. (2.1^a). In similar manner, we find 

d 


P^^^ (y I = 


d(a^ ) 


p*'^^ (y|a^) 


p(y Ic^f^) 


2 \ 2 


Map 


(ap^ 


- 6^+3 

'^y 


(A. 4) 
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and 


d(ap 

8(0^)=' (o^)^ (a^)^ 

y L y y 

+ 45 ^ - 15 

y 

In agreement with Eqs. (2.15b) and (2.l6b). 


(A. 5) 


l66 

““ lim I ■■■■IIIIIIBI nil III! !■■■ | 



APPENDIX B 


DERIVATION OF EQUATION (2.23) 

Equation (2.23) Is derived by substituting Eq. (2.20) 
Into Eq. (2.19c) and performing the resulting Integration: 


P(y) = 


1 2a 

e 


^ ^ 


-ya^/a^ 

X e y y da^ 


2 y-1 

y (yaVcjp 

2r(y) y y 


(B.l) 


Let us define 


^ ^ -X- 

y 


y 


hence , 


(B.2) 


(35 = do^ 


(B.3) 


Substituting Eqs. (B.2) and (B.3) Into Eq. (B.l) yields 

2 


y 


p(y) = 


1 


V 2 


2(^/y)C 


7T ( o'y/y ) r ( y ) ® 


Let us now define the normalized variable 


(B.it) 


n = ^ 


(B.5) 
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since d.n 



the probability density of r\ is 


P(b) 


/2tt/y r ( y ) 


oo ^ 


Yn^ 

25 5 T- 3/2 ^-5 


/2tt/y T(y) 


Y 


if 2 


r) 


Y-1/2 


/2 


K 


y- 1/2 <^’^2Y[n[) 


/2y/tt , /—I , yY-1/2 

(/2y h )' 


2^ ^^^r(Y) 


\~1/2 ^>^2Y|nl) 


(B.6) 


where Eq. (3.471.9) on p. 3^0 of Gradshteyn and 
was used in going to the second line, and where 
the modified Bessel function of the second kind 
Equation (B.6) is the same as Eq. (2.23) in the 


Ryzhlk I 50^ 
Kn ( • ) is 
of order n. 
main text. 
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APPENDIX C 


CLOSED FORM EXPRESSIONS FOR Kj^(x) AND p^(n) IN 
TERMS OF ELEMENTARY FUNCTIONS 

Prom Eq. (4) on p. 108 of Relton [5i], we have the 
recursion relation 

K ,-,(x) = — K (x) + K .(x) , (C.l) 

n+1 X n n-1 ’ 

whereas, from Eq. (8) on p. 109 of Ref. 47, we have 



Combining Eqs. (C.l) and (C.2) gives 
K3^^(x) = + i) Kj^^(x) . 

Combining Eqs. (C.2) and (C.3) with Eq . (C.l) gives 


(C.2) 


(C.3) 


= (l + ^ + ^,-1 K,Jx) 


X 


(C.4) 


In like manner, combining each previous two values of Kn(x) 
with Eq . (C.l) yields, successively: 


■'*>*> ■ k i * p • p) 

■ (l * ¥ ♦ “ * ^ ^ 

\ x x x 

Ku/a (l + ^ + M 5 + 120 ^ Ml ^ 945] 

\ ^ x2 x3 X- xW ^ 


K, (x) 


(C.5) 

(C.6) 

(C.7) 


(x)= 1 + 


^ + 210 ^ 1260 ^ 4725 ^ 10395 + 10395 I 


X 


X 


X 


X 


7 



K 


15/2 


(x) = /l + — + 2 Ii + 3150 + 17325 ^ 62370 

\ ^ X** X® 


+ 135135 + 135135 


(C.9) 


K (x) = (1 + Si + ^ + 6930 + 51975 + 270270 

\ X „2 „3 „4 „5 


X 


X 


+ 9^5945 + 2027025 2027025 \ g 

X® x’ X® / ^ 


(C.IO) 


K 


W 2 


'"l + 15 + 990 ^ 13860 ^ 135135 + 9459^5 

1 X 2 5 a f 5 


X ■ 


X 


X^ 


X ' 


4729725 ^ 16216200 ^ 34459425 ^ 34459425 \ 


X 


X 


X 


(C.ll) 


where Kn (x) is given in terms of elementary functions by Eq. 
(C.2). " 

Combining Eqs. (C.2), (C.3), (C.5), (C.9), and (C.IO), 
successively, with Eq. (2.23), we obtain the following ex- 
pressions for Py(n) defined by Eq. (2.23):' 


P, ( 0 ) = — 
' /2 


(C.12) 


Pj ( 0 ) = ^ (2 1 n i+1) e ^ I 


(C.13) 
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p (n) = — 7^ — (x^ + 6x^ + 15x + 15) e ^ , 

^ 2’*x3 

where x = v'Flnl (C.l4) 

PjCp) = ■^22560 * ^TSx'^ + 3150X- 

+ 17325X® + 62370x^ + 135135X + 135135) e“^ , 

where x = /T^|n| (C.15) 

pjn) = — (x® + 36 x^ + 630x^ + 6930 x^ 

+ 51975x'^ + 270270x^ + 945945x^ 

+ 2027025X + 2027025) , 

where x = /r8"|ril and 8! = 40320. (C.l6) 
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APPENDIX D 


ALTERNATIVE SERIES REPRESENTATION OF EXPANSION 
COEFFICIENT OF THE NONGAUSSIAN TERM 


Writing V 
( 2 ) 

efficient ]j ^2 ' 

y 


for f in Eq. ( 2 . 58)5 
in Eq . (2.24) as 


we may express the co- 


( 2 ) 

y 


u 2 = 


<J>^2(v) ly, ^(v) 


f 


h, z 


dv 


Let us define 


(D.l) 






r ^ 

^h,z 

^00 




(D.2) 


which is the "autocorrelation function" of the deterministic 
system characterization y^^ 2 ^^^ defined by Eq. (2.54). Then, 

from Wiener’s theorem — e.g., p. 54 of Ref. 34 — it follows 
that 


= I gX2TTVt 


(D.3) 


Furthermore, from Eqs. (2.8l), and (D.3), it follows by 
applying the generalized form of Parseval's theorem to Eq. 

(D.l) that we may express t ^2 as 


(2) _ 

.2 

y 


t^2 = 


cf)^ 2 (t) (f) (t) dt 

^ f T 


= 2 I <()^ 2 (t) <t>^(t) dt , 


(D.'O 
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where the second line follows from the fact that 4) 2 (t) and 

^f 

(py(t) both are necessarily even functions of t. 

Let us now consider situations where a^^Ct) fluctuates 
little over durations comparable with the duration of the 
aircraft Impulse response h(t). In these situations, the 
nominal correlation time of a|>(t) is large in comparison 
with that of z(t) — as is illustrated in Fig. D.l. Ex- 
amination of Pig. D.l suggests that we represent (p^z(t) by a 

Of 

low-order polynomial over the range 0<^T<_T , where Ty is the 
time interval over which 4)y(t) is not negligible for t>0: 


n 

<l>„2(t) = I b.f 
f j = 0 J 


0<t<T 


■y 


(D.5) 


where the bj may be interpreted as the one-sided Maclaurln 
expansion coefficients of (p zCt), 


1V(0+) 


b. = 
J 


(D.6) 




where (p 


1V(0+) 


denotes the Jth "right-hand” derivative of 


f 


4) 2 (t) evaluated at the origin. Substitution of Eq. (D.5) 
f 

into Eq . (D.4) yields the desired series representation of 

( 2 ) 

= 

y 


n 


= 2 I b, 


y 


i-o J 


t*^ 4>y(t) dt 


(D.7) 


Using Eq . (D.6), we see that the first term in the ex- 
pansion, Eq. (D.7) 5 is 
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FIG. D. 


1. BEHAVIOR OF (p 2(t) AND (p^(t) NEAR t = 0 FOR CASE 

u ^ Y 

WHERE VARIATION IN aj(t) IS SMALL OVER DURATION 
OF h(t). ^ 
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2 b 




dt = 


<P^2(0) 




dt 


Cl)„2(0) 

f 




dt 



(D.8) 


where the second line was obtained using Eq. (D.2), and the 
third line was obtained using Eqs. (2.87) and (2.92). Com- 
paring Eqs. ( 2 . 96 ) and (D.8), we see that the first terms in 
the series of Eqs. ( 2 . 96 ) and (D.7) are the same. However, 
the remaining terms differ because of the presence of odd 
powers of j in Eq . (D.7). 

Although the first correction term to the term repre- 
sented by Eq . (D.8) in the series of Eq . (2.85) had a very 

satisfying interpretation as we showed in Eq . (2. 99) 3 the 
series, Eq. (D.7), generally will be better behaved in 
practice because it utilizes the "one-sided" expansion of 
4 ^^ 2 (t) shown in Eqs. (D.5) and (D. 6 ) which does not require 
f 

that 


of E 



dt 

whereas, it is 
see Eq. (2.93) 


continuous at the origin. Thus, the existence 

is not required by the expansion of Eq. (D.7), 
required in the expansion of Eq. (2.85) — 


Finally, we note that the "one-sided moments" of (j>y(t) 
in Eq. (D.7) can be evaluated from the derivatives of the 
unilateral Laplace transform of (j)y(t) if that transform 
can be evaluated in closed form. Also notice that the 
characterization of 4)^2 (t) used in Eq. (D.5) is its one-sided 

Of 

power series expansion indicated at the very end of Sec. 1. 
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APPENDIX E 


PROOF THAT REAL AND IMAGINARY PARTS OF FOURIER SERIES 
COEFFICIENTS OF A PERIODIC RANDOM PROCESS ARE 
UNCORRELATED FOR m?«n, WHERE m,n>0 

In Eq. ( 3 . 5)5 we have expressed the complex Fourier 
series coefficients of a periodic random process 


w(t+pT) = w(t) , p=0,±l,±2^ 


by 


1 

^n T 


T/2 


•T/2 


N ~12Trnt/T 
w(t) e dt 


= a - lb_ , 
n n ’ 


where a and h are real. Consider 
n n 


(E.l) 

(E. 2a) 
(E.2b) 


c 


c* 
n m 


(a -lb ) (a +lb 
n n m 


m 


) 


and 


(a a +b b ) + l(a b -a b ) 
n m n m n m m n 


c c 
n m 


(a a -b b ) - l(a b +a b ) 
nm nm nm mn 


(E.3) 


(E.4) 


Taking the expected values of Eqs. (E.3) and (E.^), we find 
that If 


E{c c*} = 0 and E{c_c } = 0 


n m 

then we must have 


n m' 


E{a a. } ~ -E{b b } = E{b b } = 0 
n m n m n m 


(E.5) 


(E.6) 


and 
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(E.7) 


E{a^b^} = E{a^b^} = -E{a^b^} e 0 . 

Satisfaction of the two conditions in Eq. (E.5) therefore 
guarantees that all pairs of the real and imaginary parts of 
the complex coefficients c^ are uncorrelated. Cf. Davenport 
and Root pp . 91j92. 

From Eq. (E.Sa)^ we have 


E{c c*} = — 
n m ^2 } 


T/2 rT/2 


E{ w(t j ) w(t 2 )} e 


-12tt (nt j-mtz )/T 


dt jdt 2 


-T/2 -T/2 


T2 


T/2 fT/2 -12tt (nt -mt „ )/T 

J e 

-T/2 -T/2 


dt^dt^ . 


(E.8) 


If, in the inner integral, we transform t^ to x using 


T = t2 - tj ; 

hence, = t + t^,we have 


(E.9) 


-1277 (nt j -mt 2 )/T -1277 ( n-m ) t ^ /T i277mT/T 

e = e e 

therefore , 


(E.IO) 


E{c c*} = 
n m 


1. -1277 (n-m)tyT , s i277mT/T 


m2 J 

J- _t/2 


J W 

-T/2-t 


)e"'^ " dxdtj 


^ fT/2 -1277 (n-m) t J /T 

T J ^ 

-T/2 


dt J 


T 


T/2 


-T/2 


. / N 1277 mt/T , 

4>^(t) e dx , 


(E.ll) 
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where the second line is a consequence of the fact that the 
autocorrelation function of a periodic random process 

is Itself periodic with the same period T — l.e., Davenport 
and Root [35], p. 91 — and the fact that the inner integral 
in the first line in Eq. (E.ll) is Independent of ti, since 
its integrand is periodic with period T. Furthermore, for 
m 7 ^ n, the first Integral in Eq. (E.ll) is identically zero; 
hence , 


E{c c*} = 0 , m¥-r\ , (E.12) 

n m ^ ’ 

which is the first of our two conditions in Eq . (E.5). 

To check the second condition, we note from Eq. (E.2a) 
that Cjn = hence, from Eq. (E.ll), we have 


E{c c } = 
n m 


E{c c* } 


n -m' 


1 

T 


fT/2 -i2TT (n+m) t ^ /T 

] " 

-T/2 


X 



-T/2 


(J) ( T ) 

W 


^-ilrmT/T 


dx 


Consequently, 


(E.13) 


E{c c } = 0 , m/-n . (E.l4) 

n m ’ 

It follows that Eqs. (E.6) and (E.7) are satisfied for all 
m,n^0, provided m / n, which is what we sought to prove. 

As pointed out in the main text, when w(t) is generated 
from a stationary Gaussian process, the entire set of 
variates is jointly Gaussian since Eq . (E.2a) is a 

linear transformation of w(t). Hence, provided m / n, for 
all m,n>0, all aj^’s and b^^’s are statistically Independent. 
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APPENDIX F 


METHOD FOR SOLUTION OF EQUATION (3.26) FOR 
INTEGRAL SCALE OF von KARMAN TURBULENCE 


The likelihood equation for the integral scale of tur- 
bulence with negligible low frequency component Ws(t) Is 
given by Eq. (3-26): 


N 


I 

1=1 


d 

dL 


InF. (L) 
1 


r S. , N S . I 

1 _ i Y 

F. (L) N F. (L) 

- J J _ 


0 . 


(F.l) 


where we have reversed the roles of 1 and j In comparison with 
Eq. (3-26). Following the notation of Eq. (3-30), we define 


G^(L) 


^ A. 

dL 


jlnF^(L) 


(F.2) 


Let us now substitute Eq. (F.2) Into Eq . (F.l) and define 

after minor rearrangement 


E(L) 


A 1 Y 
1 = 1 


G^(L) 


N 

I 


P. (L) 


F.(L) 


(F.3) 


Then^ according to Eq. (P.3), the value of L that satisfies 
Eq. (F.l) is the value for which E(L) = 0. 


To illustrate a method for obtaining the solution 
E(L) = 0 of Eq . (F.3), we consider a vertical turbulence 
velocity record that is assumed to obey the von Karman 
(transverse) power spectral form. Prom Eqs . (3-3^) and (3-35) 
we have for the von Karman transverse spectrum: 

1+188. 75L^k^ 

F^(L) = — - ^ (F.^) 

n + 70. 78L^k^ ] 


and 
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(F.5) 


, 117.97L^k?(l-l88.75L^k?) 

G. (L) = — - ^ 

^ ^ (l+70.78L^kp(l+l88.75L2kp 

The solution E(L) = 0 of Eq . (F.3) is most easily obtained 
by trial and error. In carrying out the solution for the 
vertical record under conslderat ion ^ we used a value for the 
uppermost wavenumber corresponding to j = N in Eq. (F.3) oF 
kN = 3*0 X 10“^ cycles/meter which yielded a value of N = 6326 
points in the summation in Eq . (F.3)* 

The actual value of L = 309*^ m was obtained as follows. 
First, a trial value of L = 305 m (1000 ft) was chosen. The 
value of E(L) for this trial value of L was then computed and 
stored using Eq . (F.3)* Since the value of E(L) obtained was 
negative, a second (larger) trial value of L was chosen which 
was 335*5 meters (1100 ft). Using Eq. (F.S), a new value of 
E(L) was computed and stored using this second choice of L. 

We then possessed two values of E(L) corresponding to the two 
trial values of L. Linear Interpolation then was used to esti- 
mate a new (3rd) value of L corresponding to the value of 
E(L) = 0. For this third value of L, the true value of E(L) 
was then computed using Eq . (F.3)* The resulting value of 

E(L) had a positive sign; hence, a new (fourth) trial value 
of L was chosen which was 3*05 m (10 ft) smaller than the 
third trial value of L. The value of E(L) corresponding to 
this fourth trial value of L was then computed and stored 
using Eq. (F.3)* Finally, from the third and fourth trial 
values of L and the corresponding values of E(L), a fifth 
value of L was computed by linear interpolation corresponding 
to the value E(L) = 0. This fifth value of L = 309*^ m was 
used as the solution to the likelihood equation (F.l). These 
five values of L and the corresponding values of E(L) are 
plotted in Pig. F.l, which shows the local very nearly linear 
behavior of E(L) as a function of L. 


180 



I8I 






APPENDIX G 


TRADE-OFFS BETWEEN CHOICES OF AND m 

Here, we describe considerations to be taken into account 
in choosing values of Ch ^ the constrained least- 
squares estimation procedure described in Sec. 4. Upper 
limit in Eq. (4.5) determines the interval 0 ^ £ Ch over 

which the parameters a^, L, and a^ to a^^ in the autocorrela- 
tion function model of^Eq. (4.1) are obtained by minimization 
of the Integral squared error E in the constrained least- 
squares fit procedure. The parameter m is the degree of the 
polynomial in Eq. (4.1) that is used to represent the auto- 
correlation function of the low-frequency turbulence component 
Wg(t) over the Interval 0 <_ ^ £ 

Intelligent choices for values of Ch ^ to te used in 

the minimization procedure are not generally independent. One 
reason for this lack of Independence is the fact that our 
re present at ion in Eqs . (4.1) and (4.5) of the von Karman 
component (j)p^(^;L) of the autocorrelation function is not 

orthogonal with our representation of the low-frequency com- 
m 

V* 1 

ponent I a.^ of the autocorrelation function over the 
1 = 0 ^ 

Interval 0 £ 5 £ Ch* Thus, the low-frequency component auto- 

m 

correlation function representation a.^. has the potential 

i = 0 ^ ^ 

for representing a portion of the von Karman component of 
the empirical autocorrelation function R(5) in Eq . (4.5) in 
the Integral squares sense. However, if for given values of 
m 

3-nd m, ^ a.^. can represent, exactly, the low-frequency 
1 = 0 ^ ^ 

component R(C) over 0 £ ^ £ and if the "fast" component 

of R(^) has exactly the appropriate von Karman form for some 
values of L and ap, then this lack of orthogonality between 
m 

and \ a.^. will not be a problem. Thus ^ our goal 
1=0 ^ ^ m 

should be to choose values of and m so that \ a. 5- oan 

1=0 ^ ^ 

do a good job of representing the low-frequency component of 
the empirical autocorrelation function R(^) , while simultaneously 

m 

attempting to minimize the capability of \ a.^. to represent 

i = 0 ^ ^ 

the von Karman component of R(^). 
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To aid in implementing the above italicized rule, let us 

m 

now consider the capability of a.^^ to represent the 

i=0 ^ 

von Karman component cr^4)j^(^;L) of the autocorrelation function 
over the interval 0 £ ^ Thus, for both the von Karman 

transverse and longitudinal autocorrelation functions cl>K(^iL), 
we determine the set of coefficients aj_, i = 0,1, ••*,m that 
minimize the integral-square "error" 


E ^ 




m 


[a^<t.j^(5;L) - I 


1 = 0 


dC 


(G.l) 


To find the a.j_ that minimize E, we differentiate Eq. (G.l) with 
respect to aj , j = 0,1, •••,m which yields 


9E 

3a^. 


= -2 


m . . 

[a^4)j^(5;L) - I a.5^]C^dC. 


(G.2) 


1=0 


The solution to the set of equations (9E/3aj^) = 0, J = 0,l,*-*,m 
determines the set of aj that minimize E. From Eq, (G.2), this 
set of eouatlons can be'^wrltten as 




m 

1 


i=0 


“5J(i)j.(C;L)di; - I af = 0, j = 0,1, (G.3) 


or 


i+j+1 


’H 


m 

1=0 


a . = 


- ,^2 


Ch 


a- 1 (p^(^;L)dC, J 


= 0,1, •• • ,m, (G. 4 ) 


which is a set of m+1 linear simultaneous algebraic equations 
for the aj_ , i =0,l,«**.m in terms of von Karman auto- 

correlation function a^(p}[(^;L). These a^ minimize E. 

To write these equations in normalized form, let us define, 
as before. 
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9 


9 


(G.5a,b) 




and 

4)j^(C;L) 5 , 


(G.6) 


as In Eqs. (4.l8) and (^.20). Using Eqs . (G.5) and (G.6), 
we have 5 after dividing Eq. (G.4) by 


m 


i = 0 


j i+j+1 

2H 

i+j +1 



_i_ _ _ 

S‘'4.j^(5)d5, j 


= 0,1 


,m, (G.7) 


which is a set of equations for the normalized coefficients 
Liaj_/a^, For a given value of m, the normalized solutions ;to 
these equations depend on only one dimensionless parameter 
Thus, for either the von Karman transverse or longitudinal 
autocorrelation functions, we can solve the set of equations 
_^G.7) for Liaj_/a^, 1 = 0,1, •••,m for any family of choices of 


We shall find it convenient to have an explicit formula 
for the least integral-square ’'error” E in terms of the nor- 
malized solution vector L^aq/a^, 1 = 0,1, ••*,m to the set 
of equations (G.7)* Squaring the Integrand in Eq . (G.l) and 
rearranging terms, we have 


E = o' 




m 

(C;L)d^ -2a" a.j C ; L)dC + J d^. 


i 


r^H/ m 


i=0 


(G.8) 


or dividing by a'^L and introducing the notation of Eq . (G.6), 

we have 


I8i| 



F ^ L^a. f^u ^1 _ 

-- = 4.^(5/L)d?/L - 2 I 3^ I ^ $„(5/L)d5/L 

I ._^_2 T-1J\ 


a‘L { 


1=0 i 


Cp/ m L^a .1 

+ 1(1 V) <35/L 

U=0 L+ / 


;5h_, ™ L a f5 _ 

— “5 *K^5)dC 


<(.|(5)d5 -2 I 


1 = 0 a 


■5„/ m L a. .\ 

+ I “I i c ) dc , 

1=0 o=* ' 


(G.9) 


where we have introduced the notation of Eq . (G.5a) and (G.5h) 

However, expanding the last term in Eq. (G.9) yields 


/•Ch / rn L^a. mm L^a. L'^a. 

' ^ 1 1 „ r r 1 


1 = 0 Q' 


cm di = I I 


1 = 0 J=0 a O' 


1 I H ^ 1 +j^^ 


m m L^a. L'^" a . 

= I I — ^ ^ 


1=0 j=0 a 


2 i+j+1 


m Ea. m L+a. 

_ Y Y ^ _ii 

j=0 o^ 1=0 o^ 1+J+l 


= I 




J=0 a 




I? L a. , 


1 = 0 O' 




(G.IO) 
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where the second to the last line follows from Eq. (G.7); 
therefore the last two lines are valid only for the dimen- 
sionless solution vector Liaj_/a^ ^ i =0,1, •••^m that minimizes 
E. Substitution of Eq. (G.IO) into Eq . (G.9) yields 


E 




m 


I 

1=0 



f 1 — — — 


(G.ll) 


which is valid only for the set of dimensionless coefficients 
Liaj_/a^ , i = 0,1, •••■.m that is the solution to Eq . (G.7) — 
i.e., the set that minimizes E. 

If all = Ej we have from Eq. (G.l), 


E 

a^L 


()>^(C;L)cig/L 




(G.12) 


where the right-hand equality follows from_Eqs, (G.5) and (G.6). 
Hence, let us define a normalized "error*' E as 


E/a^ L 


(G.13) 


Furthermore, let us define a normalized set of coefficients by 



(G.14) 


Using this latter definition, our set of equations (G.7) for 

the a. 's can be written as 
1 

_ _ 

J j = 0,1, , (G.15) 

0 


m 

I 


7 i+J+1 


• ^ n i+J t- 
i = 0 - 'J - 
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whereas dividing Eq. (G.ll) by 


r%- 

(|)^(^)d^ and using the 


definitions of Eqs.__(G.13) and (G.l^), our equation for the 
normalized "error" E becomes 


m 


I a 




1 = 0 




E = 1 - 


— 


CG,l6) 


For either the von Karman transverse or longitudinal 
autocorrelation functions Eqs. (G.15) and (G.l6) 

determine the normalized "error" E In the least -squares best 
ra 

fit of I a.C to the von Karman autocorrelation function 
_ _ 1=0 ^ _ _ 

This value of E Is dependent on cholc^es of and m. 
Figures G.l and G.3 show this dependence of E on and m for 
the von Karman transverse aut ocorrelat ion function, and 
Figs. G.2 and 0.4 show the same dependence for the von Karman 
longitudinal autocorrelation function. _ Aaoordi-ng to our above 
-italicized stateme-nty large values of E are desirable . There- 
fore, for several representations 

m 

4'(5) = a^cj>„(5;L) + I a.p , 0 < C < (G.17) 

I n 1 n 


of an empirical autocorrelation function R(^)y all having 
appro ximatel-y the same capability for representing the low- 
frequency component of R(^)y but differing in values of fjj 
and my the representation with the pair of values = L^y 
and m yielding the largest value of E as determined by 
Figs. G.l to G.4 should yield the most reliable value of L. 

This rule of thumb suggests the following procedure for 
estimating autocorrelation parameters by the methodology of 
Sec . 4 . 
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1. Compute the empirical autocorrelation function R(^) — 
e.g., by the method outlined In Reference l8. 

2. By visual inspection, choose the largest value of 

m ^ 

for which \ a.^ can be expected to provide a good 

1=0 ^ 

representation of the low-frequency component of R(C) 
for each of several values of m, say m = 1,2,35 and 4. 

A different value of will generally be chosen for 
each different value of m. 

3- For each pair of values of compute L, cf-, 

a^ , a ^ , • • • , aj^ by the method described in Sec. 4. Plot 
the resulting representation, Eq. G.iPj against R(^) 
to Insure adequacy of the fit for each such computation, 

4. For each such fit — l.e., for each_palr of values of 

Ch rn — determine the value of E from Figs. G.l or 

G.3, or from Figs. G.2 or G.4, as appropriate for the 
von Karman transverse or longitudinal cases. 

5. The most reliable fit, l.e., the most reliable value 

of L, should be that corresponding to the largest 
value of E. Values of = L^h m yielding values 

of F less than, say, . 0.5 may be particularly unreliable 

m 

V 1 

because in this range l a.C has too much capability 

1=0 ^ 

for representing a portion of the von Karman component 
of the empirical autocorrelation function R(^). 
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FIG. G.l. NORMALIZED "ERROR" EQ . (G.16) IN L EAST- SQUARES 

m 

FIT OF I a. 5 '’ TO von KARMAN TRANSVERSE AUTO- 
i=0 ^ 

CORRELATION FUNCTION. 


II 




NORMALIZED "ERROR" EQ. (G.16) IN LEAST-SQUARES 
rn 

FIT OF I TO von KARMAN LONGITUDINAL 

i=0 ^ 

AUTOCORRELATION FUNCTION. 






FIG. G.3. 


NORMALIZED "ERROR" EQ. (G.16) IN LEAST-SQUAR 
m . 

FIT OF I TO von KARMAN TRANSVERSE AUTO 

i = 0 ^ 

CORRELATION FUNCTION. 


ES 
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m . 

FIT OF I TO von KARMAN LONGITUDINAL 

i = 0 ^ 

AUTOCORRELATION FUNCTION. 
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APPENDIX H 


METHOD OF APPROXIMATION OF THE INTEGRAL SCALE AND 
POWER SPECTRUM OF THE "SLOW" TURBULENCE COMPONENT w^(t) 


ExtTapolat'Lon of autoQO’PTe'lat'ion funotion model. Here, 
we develop a simple method of extrapolating the autocorrelation 
function approximation (^) of the slow turbulence component 

Ws(t) given by Eq. (4.2), 
m 

$ (S) = I a, 5-^ , 0 1 5 < (H.la) 

s j = 0 

= a +a ^+- • -+a , (H.lb) 

where the extrapolation completes the description of (^) 

s 

over the entire interval 0 £ ^ We then Integrate the 

resulting autocorrelation model to yield an approximation to 
the integral scale of the slow component WsCt); we then 
Fourier transform the resulting model to yield an analytical 
approximation to the power spectrum of the resulting model. 

The extrapolation is carried out using the simple expo- 
nential decay model 

U) = Ae““^ , 5 h - 5 

s 

which completes the range of ^ not covered by the model of 

Eq . (H.l). Parameters A and a in Eq. (H.2) are set by requiring 

(j) (F) and its first derivative to be continuous at the point 

~w 

s 

of Intersection ^ of the two models. Using this method, 

the values , • • • ,a^, and completely describe the auto- 

correlation model of Eqs . (H.l) and (H.2). 
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I 


Differentiating Eq. (H.l) gives 


0’(5)=a +2aC +*'*+ma^ 
^ 1 2 m 

s 




m 


= I Ja.S-’ ^ , 0 < Z < i 

J=1 


Therefore, we have 


m 


■J-1 


s j =1 ^ 


where the prime denotes differentiation. 
(H.2) gives 

(5) = -aAe““^H, Sjj 1 5 1 


hence, we also have 


'*'W n 


Dividing Eq . (H.5) hy Eq . (H.2) gives 


s 


Ch i 5 i 


Hence, if (^) is continuous at ^ 

■"Wg ^ 

at ^ requires from Eqs. (H.l), (H.3 


H ■ (H-3) 

(H.4) 

Dif ferentiating Eq . 

(H.5) 

(H.6) 

(H.7) 

continuity of ({) ’ (C) 

, and (H.7) that 
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m 


a = - 


• -i 

j=l ^ 
m 

y a . 

L n* 

j=0 ^ ^ 


(H.8) 


which Is an equation for a in terms of the parameters of the 
model of Eq. (H.l). To obtain an equation for A, we note from 
Eq. (H.2) that we must have 


® 

s 


(H.9) 


or. 


A = e'^^H 


m 


I 

j=0 




(H.IO) 


where Eq. (H.IO) follows from the continuity requirement of 
(^) at ^ and Eq. (H.l). Equations (H.8) and (H.IO) 

s 

yield a and A from the parameters of the model, Eq. (H.l). 

Expression for integral scale of slow component , To 
obtain an equation for the integral scale of the slow com- 


ponent, we require the integral 
and (H.2), we have 


-w 


(^)d^. From Eqs . (H.l) 


J c* 


0 

o 

11 


m 

I 

J=0 

C3 

^ J +1 

d • 

j 

J+1 

m 

y 


rJ+1 

^H 

L 

j=0 

ct « 

J 

j+1 




’H 


'H 


-a^ 


+ A 


-a 


’H 


+ A e 


(H.ll) 
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The integral scale L Is defined as 


L i 1 


(]) ( 0 


y r 


(H.12) 


Prom Eq. (H.ll), we see that an expression for A/a is required, 
which from Eqs. (H.8) and (H.IO) is 


A 

a 


\J=0 ^ / 

m . 

I jaj?/ 

j=l 


(H.13) 


Furthermore, from Eq. (H.l), we see that (0) = a„ . Hence, 

W 0 

s 

applying the definition Eq. (H.12) to (j) (?), we see from 

~ w 

s 

Eqs. (H.ll) and (H.13) that 


L = ^ 
s 0 


/ m 

m a . . ^H 


J=0 


j+1 


ra 


I ja-4"^ 

j=l 


(H.14) 


which is an expression for the integral scale L of the slow 

w 

s 

turbulence component Wg(t) in terms of the autocorrelation 
function parameters computed by the method described in Sec. 4 

Expression for power spectral density of slow component. 
By forming the Fourier transform of the above extrapolated 
autocorrelation function model, we obtain an expression for 
the power spectral density of the slow component of turbulence: 


pOO 

$ (k) = <p ( ^)cos (27rk^)d^ 

w ~w 

s s 


= 2 


4) ( ^ ) cos ( 27Tk^ ) dC 


(H.15) 
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Introducing Eqs. (H.l) and (H.2) into Eq. (H.15) and carrying 
out the integration, we have 


^ r?H ■ r 

^ (k) = 2 a . C'^cos(27rk^)d?+2A e'^'^^cos ( 2Trk^ ) d^ . (H.l6) 

j=0 J 


■H 


Using formula 2.633-2 on p. l84 of Ref. 50 for the first of 
the above two integrals and a result in the well known 
Burington tables for the second Integral, we obtain 


m 


0^) 

w 

s 


= 2 I a, 


j=0 J 


I i 


■ 3 -^ 

’H 


£=0 


( 2TTk ) 


sin(2TTk^^ + j £tt) 


-J ! (j) 

' ( 2 TTk )'^ 


, 1 . \ 

sin l 2 


-a6 


H 


-2A{ — — [27Tksin ( 27Tk^^ ) -acos ( 27rk^„ 3/ 

(an(2iTk)^ “ “ ) 


(H.17) 


However , 


£ ! I = £ ! ^ 




j ! 




(H.18) 


and 


i, 


= i I 


(H.19) 


Furthermore, for any sum, we have 
m .1 mm 


m J mm 

I I = I I 

j =0 £=0 £=0 j =£ 


(H. 20) 


Incorporating Eqs. (H.l8) through (H.20) into Eq . (H.17) and 

slightly rearranging the result yields the desired expression 

for <l> (k) : 

w 

s 
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i 


$ , (k) 
w 

s 


m 

2 I 

Ii2, = 0 


L j =^ 


m 

I a 3 I ^j+1 
TI-£) ! 


sln(2Tr^j^k + ^ JItt) 


(27TC^k) 


£+1 


-2 


m sln(^ jit) 

I a J! ^2/ 

Lj=0 (2iTk)'^ 


2Ae 


-a? 


H 


a^ + (27Tk) 


[27Tksln(27TCu^)“0^cos(2TrCTjl<^) ] j 

2 n n 


(H.21) 


where expressions for a and A are given by Eqs. (H.8) and (H.IO). 
Equation (H.21) is a closed form expression for the power 
spectral density of the slow turbulence component WgCt) in 
terms of the autocorrelation function parameters determined 
by the method described in Sec. 4. 
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APPENDIX I 


EVALUATION OF AN INTEGRAL 


We evaluate below the integral defined by Eq . (5.84): 


Y 


(3) 



' [l+(2 b+c)gM^ 

(l+bC")2(l+cC^)2 


(I.l) 


PP. 

Eq. 


( 3) 

To evaluate y j 
584 to 587 of Ref. 
(119) on p. 584 of 


we shall use the method described on 
52 . Equation (I.l) is of the form of 
Ref. 52: 


fCO 

{ • • • } d^ = 27 t 1 ^ Res(a, ) , (1.2) 

•' k 

.^OO 


where the points a^ are the poles of the Integrand that are 
located in the upper half plane. The poles of the Integrand 
occur at the zeros of the denominator of Eq. (I.l), which are 
solutions to the equations 


(l+b^2)^ = 0 (1+c^^)^ = 0 . (I.3a,b) 

These solutions are, respectively, 

= +i/l/b = +1/1/C . (I.4a,b) 

Therefore, the denominator of the integrand in Eq. (I.l) can be 
expressed as 

/b' ' /c/\ /c/ 

(1.5) 

The four factors in the right-hand side of Eq. (1.5) give rise 
to poles located at = -i//b, = +l//b, = ~l//c“, 

Ck = +i//c respectively. All four of these poles are of order 
tuo. Thus, there are two poles in the upper half plane which 
are located on the imaginary axis at 5 ^ = ±//b and = i//c 
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From p. 520 of Ref. 52 , we see that for a pole of order m, 
we require the residues 


Re s ( z . ) = 


1 


j (m-1)! 


51 


z = z . 
J 


( 1 . 6 ) 


where for our application m = 2 and 

[l + ( 2b+c ) z ^ ^ 


f(z) = 


b^c^lz^ + /z^ + 


Thus 5 we have 


(1.7) 


(z-z^)^f(z) = 


[1+ ( 2b + c ) z ^ ^ 


b ^ /z + 


1 /__2 . 1 


/b 


z ^ + — 
c 


(1.8) 


and 


(z-z^ ) ^f (z) 


[1 + ( 2b+c )z ^ ] 


2 “1 2 


b ^ c ^ f z ^ 




+ 


/ c 


(1.9) 


Differentiating Eqs . (1.8) and (1.9), as required by Eq. (1.6), 

and evaluating the resulting expressions at z = z = — and 


z = z, - — respectively, we have after slmplication 
/c 


/b 


Res ( z ^ ) = ~i 


. (b+c) ( 7b ^+c ^ ) / /b 


4b 


c-b 


(I. 10) 


and 
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/ c 
c-b 


3 


( 1 . 11 ) 


Res(z„) = 1 b(b^+bc+2e") 


Therefore, from Eqs. (I.l) and (1.2), it follows that 


Y 


(3) 



2 

7Tl[Res(z^ )+Res(z^ ) ] 


77 ^ 

(b+c ) ( 7b ^+c ^ ) 

b (b ^+bc+2c ^ ) 1 

(c-b)’ r/ 

4b'"^ 

o 

1 


( 1 . 12 ) 


which is the result given by Eq . (5.85) in the main text of the 

report . 


APPENDIX J 


EVALUATION OF SECOND ORDER PARTIAL DERIVATIVES OF JOINT 
GAUSSIAN PROBABILITY DENSITY WITH RESPECT TO ITS PARAMETERS 

The Joint Gaussian probability density of Eq . (6.4) can 

be expressed as 


p(y,y|c^,^|.Pyy) 


p 


1 B 

— ^ ^ 
2ttA " 


(J.l) 


where 


A 


A 


a"a?- 

y y 


y 


yy 


and 


4 iL 

2A 


2 ( . ) 
y y yy 


(J.2) 


(J.3) 


where 


N = (J.4) 

Denoting derivatives by superscripts in parentheses, we have 
by differentiating Eq. (J.l), 

P<1> . p(b<1> - 4^) U.5> 

and 

+ I . 

(J.6) 


P 


( 2 ) 


= P 


f. 


A 


( 2 ) 


4a 
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By differentiating the left-hand portion of Eq. (J.3)a we 
have 


B 


( 1 ) 


N 

2A 


N 


( 1 ) 


A 


( 1 ) 


N 


A 


(J.7) 


and 


B 


( 2 ) 


iL 

2A 


N 


( 2 ) 


N 


NA 


+ 2 


A 


( 1 ) 


A 


( 2 ) 

1 ' 


(J.8) 


Finally, combining Eqs . (J.7) and Eq. (J.8) with Eq. (J.6), 
we obtain an expression for the second derivative of 
p(y sY I cTy ,cr| ) in terms of the derivatives of A and N taken 
with respect to the same variable: 


(2) _ p 


2A' 


{i (*‘^>) 


- A 


( 2 ) 


A+AN^^^-3A^^^N^^^ 


+ 




- A 


(2) 


N + 


N 


(D) 


A 


( 1 ) 


A 




2A‘ 


N' 


(J.9) 


Equation (J-9) is a general expression that can be used to 

, ^ (2,0,0) (0,2,0) . (0,0,2) ^ 

evaluate p , p , and p ? ^ for use in Eq. (6.10) 

by differentiating A and N with respect to the appropriate 

variables Cy, a^, or Py^ as indicated by the superscript 

notation defined by Eq. (6-11). A and N are defined by Eqs. 

(J.2) and (J.4) . 


Let us turn now to evaluation of the cross-partial deriva- 
tives required for the last three terms in Eq. (6.10). These 


( 1 , 1 , 0 ) ^( 1 , 0 , 1 ) 


( 0 , 1 , 1 ) 


terms are and p ^ ^ .In this case, 

we shall use a double superscript notation to denote partial 
derivatives with respect to whatever two variables are required 
to evaluate these terms. 
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We can Immediately write from Eq. (J.5) 


= P 




(J.IO) 


and differentiating this expression with respect to the second 
variable, we find 


= p I- i 


a(1>1) 3a(1>0)^(0,1) 

2 “ 4a 


+ I ho, 0)3(0, 


g(l,0)g(0,l)^g(l,l) [ 


(J.ll) 


From Eq. (J.7), we can immediately write 


B 


(1,0) ^ JL / 

2A \ N A 


( J.12) 


and 


B 


(0,1) _ iL N 
2A V 


( 0 , 1 ) 


N 


A 


(J.13) 


Differentiating Eq . (J.12) with respect to the second variable 

nO,l). 

gives us B : 


pi,i) i|ui,i) APObOO) 

^ " 2A I A A 


+ 


2A 


(1, 0)^(0, l)j^ 


N 


A' 


A 


(J.14) 
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Finally, substituting Eqs . (J,12) to (J.l^) into Eq. (J.ll) 
yields the desired expression for 


(1,1) ^ ^ )3 ^(1, 0)^(0, l)_^(l,l)^^^j^(l,l) _ 3 ^(0,1),, (1,0) 

2)2 2 


2A‘ 


- I + 


3a(1,0)a(0,1) ^ (1,1) 

A 


N 


^(l30)^(0,l) 


^(0,1) qN 


(0 1) A 
2k " + - 


(1,0) (0,1) 



2A^ 


(J.I 5 ) 


Equation (J.15) is a general expression that can be used to 

. , (1,1,0) (1,0,1) . (0,1,1) ^ ^ 

evaluate p , P j and p ^ ^ for use in Eq. 

(6.10). In evaluating the various terms in Eq . (J.15)j the 

double superscript notation is used to denote derivatives with 

respect to whatever two variables in and 

p(0,l,l) partial derivatives are taken with respect to. 

A and N are defined by Eqs. (J.2) and (J.^). 


As a check of Eqs. (J.9) and (J. 15 ), we note that the 

expression for should reduce to the expression for p^^^ 

when we substitute the right-hand sides of the following 

expressions into Eq . (J. 15 ) 


a(1,0) ^ ^(1) 


^(0,1) ^ ^(1) 


A^l^l) = A^^^ 




j^(0,l) _ 




Carrying out these substitutions reduces Eq . (J.15) to Eq . (J.9). 


We shall now use Eq. (J.9) to evaluate the terms 
p(0,2,0)^ and for use in Eq. (6.10), and following 

that we shall use Eq. (J.15) to evaluate 

^ ( 0 , 1 , 1 ) 

and p ’ ’ . 

To evaluate p ^ ^ ^ ( y ,y I , a|-, ty^) s we identify all 

derivatives in Eq . (J.9) as derivatives'^ with respect to 
[as indicated by Eq. (6.11)]. Therefore, from Eqs. (1.2)“^ 

and (J.4) we have for evaluation of ^ ^ ^ ^ ^ ^ (y cr?, \iy^) : 
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0 


(J.l6) 


A 


(1) = 



A 


(2) ^ 


= -y^ 


= 0 


(J.17) 


When Eqs. (J.l) to (J.4) and Eqs . (J.l6) and (J.17) are sub- 
stituted into Eq. (J.9) we obtain the desired expression for 


p(2,0,0) 


(y,y 







In like manner, we have developed the following table 
which gives the evaluates of A^^\ m(2) 


and N 


required for use in Eq. (J.9) to obtain expressions for 


(2,0,0) p(0,2,0) 


and p 


( 0 , 0 , 2 ) 




a(2) 



p(2,0,0) 

cl 

0 

-y^ 

0 

y 



p(0,2,0) 

y 

0 

-y^ 

0 

p(0,0,2) 

"^^yy 

-2 

2yy 

0 


Table J.l. Evaluates of A^^^, and for sub' 

stltutlon into Eq. (J.9) to determine expressions for 

p(2,0,0)^ p(0,2,0), p(0,0,2)^ Quantities p. A, B, and 

N are given by Eqs, (J.l) to (J.4) respectively. 


Similarly , 
and p 


we obtain 

( 0 , 1 , 1 ) 


Table J.2 for evaluation of 
Eq. (J.15): 
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p(l,l,0) 

^(1,0) 

a(0,1) 

y 

^(1,1) 

n(W0) 

„(0,1) 


a? 

y 

1 

. 2 

-y 

-y^ 

0 

pdmu) 

y 

-2d . 

yy 

0 

-f 

2yy 

0 

p(0,l,l) 


-2d . 

0 

2 

-y 

2yy 

0 


y 

yy 


1 




Table J.2. Evaluates of 

and for substitution into Eq. (J.15) to deter- 

4 > ( 1 , 1 , 0 ) ( 1 , 0 , 1 ) , ( 0 , 1 , 1 ) 

mine expressions for p , p , and p ’ ’ . 

Quantities p. A, B, and N are given by Eqs . (J.l) to (J.4) 

respectively . 

When Eqs. (J.2), (J.4), and the parameter values given in 
Table J.l are substituted into Eq . (J.9)j we obtain the 

V, ^ (2,0,0) (0,2,0) . (0,0,2) . , 

expressions for p , p , and p ’ given by 

Eqs. (6.15) to (6.17) respectively, and when Eqs. (J.2), (J.4), 
and the parameter values given in Table J.2 are substituted 

into Eq. (J.15)j we obtain the expressions for 

pdaOjl)^ and given by Eqs. (6.l8) to (6.20) respec- 

tively. Since these derivatives are shown evaluated at the 
expected values of the parameters Oy , , and bars are 

shown over these parameters in Eqs. (6.15) to (o.20); further- 
more, we have used the fact shown by Eq. (6.?4) that 

y.EE{y.}=0 (J.l8) 

yy "^yy 

in the expressions of Eqs. (6.15) to (6.20). 


207 



APPENDIX K 


DERIVATION OF INPUT-RESPONSE RELATIONS FOR 
INSTANTANEOUS CROSS-SPECTRAL DENSITIES OF 
NONSTATIONARY STOCHASTIC PROCESSES 


Here, we derive input-response relations for instantaneous 
cross-spectral densities that are direct extensions of the 
results derived in Ref. 3^. Let us define the instantaneous 
cross-correlation function of two real, generally nonstationary 

stochastic processes {x.(t)} and {x, (t)} as 

J 


h.x, = 

J L 


E { X . ( t 
J 




+ 


p}. 


(K.l) 


When the two processes are identical — l.e., when j = k, and 
therefore xj(t) = x^(t) — c{) ^ (T,t) is an even function of t 

j 

as is immediately apparent from Eq. (K.l). However, when 
x.(t) 7 ^ X, (t), (j) (T,t) is not generally an even function 

J K X . X, 

J k 

of T. The definition, Eq. (K.l), is a direct extension of 
the definition, Eq. (7), of Ref. 34. 

We define the instant aneous cross-spectral density 
^ (f,t) of the two processes {x.(t)} and {x, (t)} as the 

J k 

Fourier transform with respect to t of cf) (x,t) — i.e., 

X • X 


^ (f,t) = 

X . X, ^ 

J k 


. J k 


-127Tf T 


dx 


(K.2) 


which is a direct extension of the definition, Eq. (9a), of 
Ref. 3^. However, in the present case ^ (f^t) is not 

X » Xt 

J k 

generally real and an even function of f as it is in the 
case where xj(t) = xj^(t). Let us further define the Fourier 
transform with respect to t of 0 (f,t) as 

X • X-i 

J k 
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#00 

J k i„ J k 


(K.3) 


which Is a direct extension of Eq. (27) of Ref. 3^ • 

We may relate (f,v) directly to the Fourier trans- 

J k 

forms of the sample functions xj(t) and Xj^(t). Substituting 
Eqs. (K.l) and (K.2) Into Eq . (a. 3) and Interchanging the 

orders of expectation and Integration, we have 


»oo #oo 


J k 


x.(t- ^)x^(t+ |-)e "^^^dTdtV. 


(K.4) 


— GO — .OO 


Let us now transform to the new variables of Integration 


hence , 


t = t — — 
'^1 ^ ? 


T = t^-t^ 


t, = t + ^ ; 


t = 


tj+t^ 


(K.5) 

(K.6) 


where | 3 ( x ,t )/3 ( t ^ , t ^ ) | = 1. Substitution of Eqs. (K.5) and 
(K.6) Into Eq . (K.4) and using 


dTdt = I a ( t , t )/9 (t J ,t 2 ) I dt ^dt 2 = dt^dt^ , 

we have after minor rearrangements 


(K.7) 




x.(t,)x. (t^)e 
1 1 k 2 ' 


-12Tr[(f- ^)t^-(f+ ^)tj] 


dt dt ^ 
1 2 


127T ( f+ ~)t 

E{| x^. (t^)e ^ dt^ 


-12n(f~ 

Xj^(t2)e dt^ 


E{X*(f+ |-)} , 


(K.8) 
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where we have defined 


X.(f) ^ x.(t)e-^2wft 

J I 0 


(K.9a) 


X,^(f) = I x,^(t )e“^^’^^^dt , 




k 


(K.9b) 


and where the superscript asterisk In Eq . (K.8) denotes the 

complex conjugate. Thus, from Eq . (K.8) we see that $ (f,v) 

X • X, 

J 

can be expressed directly In terms of an expectation of the 
cross-products of the Fourier transforms of the sample func- 
tions of the two processes {x.(t)} and {x, (t)}. 

<3 

The Fourier mate to Eq. (K.2) Is 

r°° 

‘•‘x.x " '*x.x . (K.IO) 


Combining the evaluations of Eqs. (K.l) and (K.IO) at t = 0 
give s 


(|) (0,t) = E{x.(t)x(t)} 

Xj Xj^ J K 




(K.ll) 


which Is the extension of Eq. (12a) of Ref. 3^ to Instantaneous 
cross-spectral densities. For any time t. Integration of the 
(complex) Instantaneous cross-spectral density $ (f,t) 

^j^k 

over all f gives the expected value of the product x.(t)x, (t) 
at that same Instant of time t. ^ 

Cross-spectral density input-response relations . Consider 
the response {y(t)} of a linear t Ime -Invariant system to an 
Input process {x(t)}. Let h(t) denote the unit Impulse 
response of the system. For any Input sample function x(t) 
we have 
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y(t) = 


oo 


(K.12) 


X (t ) h(t-T )dT 


3 


as in Eq . (1.10b). Let H(f) denote the complex frequency response 

of the system as defined by Eq . (1.9), and let X(f) and Y(f) 
denote the Fourier transforms of x(t) and y(t) defined in the 
same manner as in Eq. (1.9). Then, it is well known [e.g., 
p. 57 of Ref. 29 ] that X(f), Y(f), and H(f) are related by the 
product 


Y(f) = X(f)H(f) . (K.I 3 ) 

Let {xj(t)} and {xj^(t)} denote two different input processes 
and let {yjXt)} and (y]^Tt)} denote the corresponding response 
processes. That is, each sample function from process fxj(t)} 
generates a response sample function by the relation 


y (t) = h (T)x.(t-T)dT 

J J J J 

—00 


(K.14) 


and each sample function from (x^(t)} also generates a comparable 
response sample function 


y 


k 


(t) = 


h ( T ) X ( t-T )dx 

JK. K, 


(K.15) 


Impulse response functions hj(t) and h^(t) are potentially dif- 
ferent. The frequency domain counterparts of these input- 
response relat ionships are 

Y (f) = X (f)H (f) (K.16) 

J J J 

and 

Y^(f) = X^(f)H^(f) . (K.I 7 ) 

From Eqs. (K.I 6 ) and (K.17)j we therefore have 
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E{Y*(f+ = E{X*(f+ |))H*(r+ |), (K.l8) 

or, applying the definition, Eq. (K.8), to input processes, 
output processes, and system complex frequency response functions, 
we have 

*y.y, = h.x, ’ (K.19) 

J k j k j k 

where no expectation operation is required in defining ^ (f,v) 

J k 

because hj(t), hj^Ct) and their transforms are here assumed to be 
deterministic . 

Equation (K.19) expresses the transformed Instantaneous 
cross-spectral density of the response processes {y-(t)} and 
{yk('t)^ cLS the product of the transformed Instantaneous cross- 
spectra of the input processes and system Impulse response 
functions. Hence, from the Fourier mate to Eq . (K.3), 


J k 


<!> (f,v e 

X . X, ^ 

, J k 


-ifirvt 


dv 


(K.20) 


and the analogous relations for the response cross-spectra system 
impulse response cross-spectra, we have by applying the con- 
volution theorem to Eq. (K.19): 


$ (f,t) = (f,u)du . (K.21) 

J k J k J k 

Equations (K.19) and (K.21) are the Instantaneous cross-spectral 
density input -response relations that are direct extensions of 
the Instantaneous auto-spectral density input-response relations, 
Eqs . (39) and (40) of Ref. 3^. 

Reduot'ton to the case of stationavy input processes. Let 
us now consider the case where (p (x,t) defined by Eq . (K.l) 

X • X T 

J k 
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is independent of t — i.e., is dependent only on the time dif- 
ference T. In this case, we see from Eq, (K.2) that (f,t) 

yi • 

J k 

also is independent of t . Denote this independence by replac- 
ing t by a vertically centered dot. It then follows immediately 
from Eq. (K.21) that $ (f^t) also is independent of t, which 

we also shall denote by a vertically centered dot. Thue^, when 
(p (Tyt) is independent of Eq, (K,21) veduoes to 

J k 


^ (f , 


) = 


^ ( f , 

J k 


^h.h, 

. J k 


(K. 22 ) 


However , 


from the counterpart of Eq . 


(K.3) applied to h . h, 

3 


1 . e 


j 


I CO 

$ , (K.23) 

CO J k 


we have by setting v = 0 in Eq . (K.23)j 


$ (f,t)dt = ? (f,0) 

, J k J k 


= (f) 

u ^ 


(K.24 ) 


where the second line follows directly from Eq. (K.8) applied 
to hjhk rather than x.-xj^. Combining Eqs. (K.22) and (K.2^) 
yields 

v = ^x X (f'.-)HHf)H (f) . (K.25) 

yjyk x^.x^ j k 

Finally, by applying Eq. (K.ll) to the response yjyk^ we have 
from Eqs. (K.25) for stationary input processes ^ 
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CO 


(K.26) 


E{yj(t)yj^(t)} = J (f,OH*(f)H^(f)df , 

_oo J ^ 

which is valid whenever the instantaneous input cross-correlation 
function, Eq. (K.l), is independent of t. 
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